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Abstract 

Optical Tomography of the ZaP 
Flow Z-Pinch Plasma 

by Adam Matthew Madson 

Chair of the Supervisory Committee: 

Associate Professor Uri Shumlak 
Department of Aeronautics & Astronautics 

The ZaP Flow Z-Pinch Experiment forms a Z-pinch with an embedded axial flow 
that exhibits stability for 500-2000 times the instabililty growth time for a static 
Z-pinch. Two linear photodiode arrays placed orthogonal to each other measure 
chord-integrated emissivity of the plasma. Deconvolving these data produce time 
resolved 1-D and 2-D emission prohles of the Z-pinch cross-section during the sta¬ 
ble quiescent period. These unhltered emission prohles show that the pinch has a 
radius of approximately 1 cm. Filtering the light for Ha emission shows a dehned 
hollow prohle 1 cm in radius. Bremsstrahlung hltered light exhibits a peaked prohle 
during periods of high intensity emissivity and a hollow prohle at lower intensities. 
Additional measurements from an Imacon camera, interferometer, and azimuthal ar¬ 
rays of magnetic probes determine the spatial correlation between plasma emissivity, 
plasma density, and plasma current. The tomographic information collected by the 
photodiode arrays increases understanding of how stabilization theory and Z-pinch 
conhguration plasmas. 
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Chapter 1 

INTRODUCTION 

With the decline of traditional fossil fuel reserves around the globe, there is great 
interest in the research and development of new “alternative” energy. A safe, clean, 
reliable energy with a nearly unlimited or even renewable source is the objective of 
current research in the energetics held. Fusion reactions, the processes where the 
combination of light atomic nuclei into heavier nuclei release substantial amounts 
of energy, offer an exciting area of research in this held. These reactions naturally 
occur inside the sun, where predominately hydrogen ions form helium nuclei through 
the proton chain process [11]. Unfortunately, this is a slow and weakly interacting 
process, and it is only the massive central pressure inside the sun that increases the 
reaction rate to a level that sustains fusion. The gravity created by the sun contains 
this reaction, providing a counter-balance to the outward force resulting from the 
pressure gradient between the ionized gas at the core and at the surface. Therefore, 
a successful fusion reaction depends on the containment of a highly dense volume 
of mostly ionized gas, called plasma. As dehned by Francis Chen, “A plasma is a 
quasineutral gas of charged and neutral particles which exhibits collective behavior 
[3].” Several different containment schemes for highly dense plasmas are currently 
under investigation around the world, including Spheromak, Tokamak, Reversed Field 
Pinch, and Z-pinch plasmas. The ZaP experiment at the University of Washington 
Aerospace and Energetics Lab attempts to contain high density Z-Pinch plasmas 
using an innovative sheared axial flow configuration. Two multi-channel photodiode 
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arrays measure optical emissions from the plasma, comparing them to other diagnostic 
systems as well as providing one-dimensional and two-dimensional deconvolved radial 
prohles. 


1.1 Introduction to the Z-Pinch 


A Z-pinch is a relatively simple plasma conhguration in terms of shape and con¬ 
tainment requirements. It is a cylindrically symmetric plasma column with an axial 
current, Jz, as shown in Figure 1.1. This axial current induces an azimuthal mag¬ 
netic held Bg which in turn conhnes the plasma radially. Without any plasma and 
resulting current, there would be no magnetic held. Therefore, the Z-pinch is entirely 
self-contained, with no externally driven magnetic helds needed for its creation. The 
simplicity of this conhguration greatly reduces the complexity of the force balance 
equations governing the equilibrium pressure balance. 

Ampere’s law and the ideal magneto-hydrodynamic (MHD) force balance equation 
are 

W X B = (1.1) 

J X B = Vp. (1.2) 


Given that Jz and Bg are the only non-zero current and magnetic held, these equations 
reduce respectively to the following [6] 


J. = -4- i^Bg) 

dr 

T R - 

■Jztjg — —r- 
dr 


Substituting 1.3 into 1.4 yields 


£ ( 

dr V 2/io/ Hor' 


( 1 . 3 ) 

( 1 . 4 ) 


( 1 . 5 ) 


The derivative terms represent the plasma particle pressure and the magnetic pressure 
respectively. The plasma pressure force is always pointing out from the radial axis 
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Figure 1.1: Z-pinch geometry 


of the Z-pinch. The next term is the magnetic pressure force, which points outwards 
when confining the plasma outside the plasma radius. The term on the right hand 
side is the magnetic tension term, created by the curvature of the magnetic held 
lines. For a discrete plasma with a uniform current prohle out to its radius, the 
plasma and magnetic pressure forces balances the magnetic tension force outside the 
plasma column. An easy way to illustrate magnetic tension term is by visualizing the 
force generated by a rubber band wrapped around an object. It is this inward force 
that provides the radial force balance needed to conhne a Z-pinch. This balance also 



4 


results in greater confinement efficiency, represented by the term f3. 

^ average plasma pressure < p > g, 

peak magnetic pressure Bg/2po 

This conhnement makes the Z-pinch a good candidate for fusion power applications, 
unfortunately, real life instabilities also make it violently unstable. 


1.2 Instabilities Inherent to a Z-pinch 

Z-pinches are unstable to both current driven kinks and pressure driven MHD insta¬ 
bilities [6]. These instabilities arise when a displacement in the plasma offers a lower 
energy state versus the existing conhguration. Any variations in the current prohle or 
the Bg held can lead to instability. The two best examples are the m = 1 and m = 0 
instability modes. 

The m = 0 mode, or sausage mode shown in Figure 1.2(a), occurs when a small 
section of the plasma column radially contracts. The reduced cross-sectional area 
increases the magnetic held at this location and therefore its magnetic tension force 
directed inward. This results in even further narrowing of the plasma at this location, 
until the plasma physically breaks apart. 

The m = 1 mode, or kink mode shown in Figure 1.2(b), occurs when the plasma 
column bends, creating a stronger local magnetic held on the inner side of the bend 
and a weaker held on the outer side. The resulting increase in magnetic tension 
creates a force in the direction of the kink, further increasing its deformation un¬ 
til containment is no longer possible [6]. Traditional methods for stabilization have 
included modifying the pressure prohle, reducing wall radius, and applying axial mag¬ 
netic helds. Unfortunately, these methods all have signihcant drawbacks that reduce 
their usefulness for fusion plasma applications. Limiting the radial pressure gradient 
for example, can stabilize the plasma against the m = 0 mode. Specihcally, if the 



(a) 


(b) 


Figure 1.2: (a) m = 0 “sausage” mode and (b) m = 1 “kink” mode with surrounding 
Be field 


pressure profile satisfies the following relationship 


4r dlnp 

2 + r/3^ dlnr 


( 1 . 7 ) 


where F is the ratio of specific heats and (3 is the ratio of plasma pressure to magnetic 
pressure, the m = 0 mode can be marginally stabilized [15, 19]. Unfortunately, this 
has no effect on stabilizing the m = 1 mode. 

Likewise, an imbedded axial magnetic field in the plasma can connteract the mag¬ 
netic tension of the m = 1 mode, bnt it also limits the maximnm axial cnrrent and 
opens np the field lines [16, 17]. As the plasma colnmn bends, the axial field stretches, 
providing a larger and larger force to connteract the instability. Limiting of the cnr¬ 
rent lowers the maximnm attainable plasma pressure, making this option nndesirable 
for fusion power applications. These drawbacks in conventional stabilization methods 
resulted in the rednction of Z-pinch research in favor of other plasma configurations. 




6 


1.3 Sheared Flow Stabilization 

Theoretical research by Shumlak and Hartman [18] and experimental research by 
Shumlak et al. [19], provides evidence that sheared velocity flow offers the exciting 
possibility of limiting instability propagation in the Z-pinch. Specifically, a radial 
shear of the axial flow, dvz/dr > O.l/cVd, where a is the plasma radius, k the wave 
vector, and Va the Alfven velocity, can satisfy the marginally stable pressure profile 
condition previously given by Equation 1.7. Alfven velocity relates to the electro¬ 
magnetic Alfven wave. An out of plane electric field produced by the wave combines 
with the Bq magnetic held to produce an E x B drift [3]. This results in the plasma 
particles moving in the same oscillating variation as the magnetic held lines, motion 
known as frozen-in flow. Perturbation propagation occurs as these variations in the 
plasma travel down the held lines and cause local instability modes. The speed of this 
propagation is given by the Alfven velocity, defined by Va = Bo/^fxop [3]. There¬ 
fore, relating the Alfven speed to radial velocity shear allows for the construction of a 
normalized stability criterion for varying radii. This criterion is shown in Figure 1.3. 


1.4 Overview of the ZaP Flow Z-pinch 

The ZaP Flow Z-Pinch Experiment at the University of Washington Aerospace and 
Energetics Research Laboratory currently investigates the effects of sheared axial 
flow on Z-pinch plasma stability. In its present design, a coaxial acceleration region 
similar to a Marshall gun in an attached plasma assembly region produces 1 meter 
long Z-pinch plasma columns. The present design of the experiment is included in 
Figure 1.4 

To form the Z-pinch, hydrogen gas is injected into the vacuum chamber between 
the inner and outer electrodes near the midplane of the acceleration region. Applying 
a negative voltage to the inner electrode causes the gas to ionize and form a current 
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Threshold Shear 



Figure 1.3: Normalized velocity shear threshold, r^/a is wall radius over plasma 
radius, v' is dv^/dr, Va is the Alfven velocity, and k = 27r/A is the wave number. 
Plasmas whose edges are close to the wall are generally stable, however for fusion 
applications the pinch radius is much smaller than the wall radius. The threshold 
stability for this case would correspond to the right side of the hgure, with a value of 
0 . 1 . 
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Acceleration Region 


Assembly Region 


1 m 


1 m 


Figure 1.4: The ZaP Flow Z-pinch Experiment 


sheet running from the outer to the inner electrode. This current induces a magnetic 
held in the azimuthal direction, and the two combine to produce a J x B Lorentz 
force in the z direction that pushes the plasma down the inner electrode and into 
the assembly region. As the current moves into this region, the increasing distance 
and slant between the inner electrode nosecone and outer electrode connection point 
produces a component of the Lorentz force pointing towards the ^-axis, causing the 
current sheet to collapse down onto axis. Once there, the plasma pressure increases 
and the Z-pinch forms with an imbedded sheared axial how. This sequence is shown 
in Figure 1.5. 

Neutral gas remaining in the acceleration region continues to ionize and feed into 
the pinch until no more remains or instabilities destroy the pinch. Experimental 
results show that the ZaP sheared how Z-pinch remains stable 500-2000 times longer 
than the instability growth time of a static Z-pinch. Other characteristics of the ZaP 
Z-pinch are included in the Table 1.1 [14]. 
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(d) (e) (f) 

Figure 1.5: Z-pinch formation, (a) Neutral gas is injected in the acceleration region 
and (b) ionized when a potential is applied between the electrodes, (c) The J x B 
Lorentz force accelerates the current sheet into the assembly region (d and e) and 
as it extends in length, the resulting inward Lorentz force pushes the plasma onto 
axis, (f) This creates a dense Z-pinch column aligned on axis with an embedded axial 
shear. 
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Table 1.1: ZaP Flow Z-pinch typical operating parameters. 


parameter 

symbol 

value 

ion velocity 

Vion 

5-10 cm/ps 

Alfven speed 

va 

12-15 cm//is 

quiescent period length 


20-40 ps 

total temperature 

Te + Ti 

150-200 eV 

electron number density 

rie 

1016-1017 

edge magnetic field 

-Sedge 

1.0-1.8 T 

Z-pinch radius 

a 

1 cm 

Z-pinch length 

L 

1 m 

peak current 

-^p max 

150-250 kA 

peak power 

Sp max 

0.2-1.5 GW 

capacitor charge voltage 

^mb 

5-9 kV 
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Chapter 2 

DIAGNOSTIC SYSTEMS AND COMPARISONS 

The ZaP experiment currently uses a variety of diagnostic systems to investigate Z- 
pinch plasma stability. These include measurements of the magnetic held and plasma 
current using axial and azimuthal arrays of magnetic probes, plasma currents using 
Rogowski probes, chord-integrated plasma density from a two chord interferometer, 
and light intensity or emissivity from CCD and ICCD spectrometers, photomultiplier 
tube, fast frame camera, and photodiode arrays. Much time and research goes into 
the secondary conhrmation of data taken by a single diagnostic device. Relating 
photodiode emissivity data to the current location measured by the magnetic probes, 
the electron density measured by the interferometer, and the emission prohle taken 
by the fast frame camera provides insight to the stabilization theory of the sheared 
how Z-pinch. 

2.1 Photodiode Array (PDA) 

Photodiodes are relatively simple devices that prove valuable in converting light inten¬ 
sity to a current. A photodiode crystal consists of a surface P-layer and a subsurface 
N-layer. Photons striking the crystal stimulate electrons in these layers out of their 
energy levels. The electron-hole pairs move in opposite directions, with the electrons 
toward the N-layer and holes toward the positive P-layer. When placed in a closed 
circuit, the voltage diherence caused by this movement creates a current that varies 
with light intensity. As the resistance of this circuit approaches zero, the current 
varies more and more linearly with the changing light level, as shown in Figure 2.1 
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Negative Bias 



Figure 2.1: Photodiode current vs. voltage curve. In a closed circuit with no resis¬ 
tance, current changes linearly with respect to light intensity. 


[10]. The depletion layer created on the anode side of the photodiode between the P- 
and N-layers slow the response time between light intensity variations and the cur¬ 
rent. Introducing a reverse bias current on the photodiode replenishes this depletion 
layer, thereby reducing signal response time. More discussion regarding the biasing 
circuit is in Section 4.1. 


2.2 Magnetic Probe Arrays and m = 1 Mode Measurements 


The magnetic probes in the ZaP experiment measure mode activity and current prop¬ 
erties of the Z-pinch by measuring the changes in the local magnetic held. Using the 
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relationship described by Faraday’s law [9] 

^E d\ = - j B da, (2.1) 

coils of wire imbedded in the outer electrode generate a voltage based upon the varying 
magnetic held. A single axial array of probes measures the current sheet distribution 
and location. Four azimuthal arrays, three in the assembly region with eight equally 
spaced probes and one in the acceleration region with only six probes measure the 
azimuthal mode activity [8]. Normalizing the m = 1 magnetic mode by the aver¬ 
age magnetic held at the axial location of the array indicates the displacement of 
the Z-pinch. This is especially important since it allows for the creation of an em¬ 
pirical boundary condition between a stable pinch, dehned as the quiescent period, 
and unstable regions, dehned as the formation and disruption phases. The MAT- 
LAB code magneticmode .m used to determine magnetic mode activity is included in 
Appendix A.3. 

2.3 Interferometry and Chord-Integrated Electron Density 

The two chord Mach-Zehnder conhguration interferometer determines electron density 
by measuring the index of refraction with two coherent laser light sources. One beam 
travels through the plasma and then recombines with a reference beam, resulting 
in constructive and destructive interference. Photodiodes record these changes, and 
since the phase of the laser light shifts proportionally to the chord-integrated electron 
density, multiple beams traveling through the plasma at a single axial location allow 
for the determination of the line-integrated density through different radial locations 
of the plasma. 

2.4 Imacon Camera 

Taking optical measurements of the plasma is simply a matter of recording the image 
viewed through a window in the vacuum chamber onto some sort of medium. High 
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Time (jlls) 

Figure 2.2: Normalized m = 1 mode displacement. Quiescence occurs when the 
displacement is below the empirical value of 0.2 as shown by the dashed black and 
green lines. This corresponds to a radial displacement in the tank of approximately 
1 cm off the z-axis. In this example, Z-pinch quiescence occurs from 37.0 — 85.4 fis. 


speed photography offers an excellent means of accomplishing this task and storing 
a great deal of visual data on a small section of film. An Imacon Model 790 camera 
takes images and stores them on Polaroid 667, 3000 speed film. Data are stored as 
a series of framed images at stepped time increments or as a streak photo along the 
length of the Polaroid him. An electrical circuit in the Imacon controls the geometric 
placement of the image on the him through the use of “shutter” and “compensating” 
capacitor plates. For a framing picture, the camera electronically places images at 
separate locations on the him. In the case of the streak photo, a narrow slit is used to 
view only a thin vertical section of the plasma. The camera then sweeps this image 
along the length of the him at a known linear rate. This creates an image of plasma 
movement in the vertical direction versus time in the horizontal direction of the him. 
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Chapter 3 

TOMOGRAPHIC THEORY 

The use of photodiodes to collect radiated data is just one of a variety of techniques 
used in the held of image reconstruction, or tomography. It is used extensively in the 
scientihc and medical helds, allowing for an internal prohle of an object through 
the use of external viewing probes. The development of the theoretical basis for 
tomographic methods is credited to Alan Cormack in the mid 1960’s [4, 5]. His 
research, and its subsequent application by Godfrey Hounsheld resulted in the creation 
of computer-assisted tomography, or the CAT scan, and a Nobel Prize in Physiology 
and Medicine in 1979. While the CAT scan method involves radiating an object 
with an external X-ray source and combining 2-D prohles into a single 3-D image, 
the tomographic reconstruction of plasma is slightly different. Since the plasma is a 
radiating body, the measuring probe can be a passive collector instead of an irradiating 
source. The possible lack of symmetry also requires the use of multiple detectors in 
various viewing axes to deconvolve the plasma prohle. Using multiple photodiodes 
in an array introduces the possibility of generating a prohle of the visible emissivity 
along the viewing chords. The ZaP how Z-pinch makes use of two photodiode arrays 
to measure chord-integrated emissivity, as well as tomographically create one and two 
dimensional radial prohles of the plasma interior. 

3.1 Plasma Bremsstrahlung and Line Radiation 

In order to gain an accurate understanding of the plasma from tomographic data, 
it is important to determine exactly what sort of light the photodiodes detect, their 
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relationship to density and energy, and how they affect the resulting current. Energy 
radiates from the Z-pinch in many different forms, but for this experiment, the two of 
greatest interest are bremsstrahlung and line radiation. Bremsstrahlung or braking 
radiation results from the interaction of ions and electrons undergoing Coulomb colli¬ 
sions. As fast moving electrons approach the relatively stationary ions, the long-range 
Coulomb field deflects them into a hyperbolic orbit. This deflection is determined from 
the point of closest approach in the absence of any electromagnetic forces, called the 
impact parameter (&), and the deflection angle {6), governed by the following equation 

[7] 


2 {b/bo 


sin0 = ^A^^ (3.1) 

1 -|- {b/bo) 

where bo is the relative impact parameter when deflection is 90°. As the electron 
orbits the ion during Coulomb collision, it accelerates, resulting in the production of 
electromagnetic radiation called bremsstrahlung. The power radiated, W, from the 
interaction between acceleration and the Coulomb force is 


W = 




(3,2) 


3 (47reo)^ 

In this equation, Zefr is the effective ion charge, e the electron charge in Coulombs, 
eo the permittivity of free space, m the ion mass, c the speed of light, and r the 
instantaneous radius of the electron from the ion. Integrating this value over the 
interaction time produces the following amount of emitted radiation per collision. 


bhrad 




(3.3) 


3 (47reo)^ m?(^vb‘^ 

where v is the electron velocity relative to the ion (normally taken as the electron 
thermal speed). The summation of this equation for all Coulomb collisions represents 
the total amount of energy emitted from the plasma from bremsstrahlung. Manipu¬ 
lating this equation yields the power radiated per unit volume of the plasma, given 
by 




br 


31/2 (dvreo)^ 7 ip/‘^P'h 


(3.4) 
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where Tg is the electron temperature, h is Planck’s constant, and n* and Ue are the ion 
and electron number densities respectively. Therefore, bremsstrahlung from a plasma 
increases linearly with temperature and density due to the higher rate in Coulomb 
collisions per unit volume. 

Line radiation on the other hand, arises from the transitions of bound electrons 
between energy levels. The power produced per unit volume relates to the combi¬ 
nation of electron density and high-Z ion density, as well as electron temperature 
[7]. The distinct energy levels of different ions result in line radiation occurring at 
specihc wavelengths. This allows for the identihcation of various impurities inside 
the vacuum chamber of the experiment that superimpose over the bremsstrahlung 
radiation. Figure 3.1 is the radiation measured on the ZaP experiment for a 5 kV ca¬ 
pacitor bank charge. The greatest line radiation contributions come from ions such as 
N-III, O-II, O-III, Ha (H-I), B-III, and C-III. With the exception of Ha at 656.3 nm, 
most other signihcant line radiation occurs between 300-500 nm. Placing high-pass or 
band-pass Liters in front of the collimator openings of the photodiode array therefore 
significantly reduces the viewed line radiation. 

3.2 Centroid Calculation 

Given a linear array of photodiodes, it is very simple to calculate the center of emis- 
sivity of the viewed radiation from the chord-integrated data. The following equations 
govern the location of the centroid in two dimensions. 

f X dArea f y dArea 

^ Area ’ ^ Area 

Since data from the photodiodes is taken at discrete points along the array axis, it 
is easiest to determine the overall centroid by summing the centroid components of 
many smaller subsections using matrix arithmetic. In this case, the centroid is the 
summation of individual chord spacing Axi times the chord position Xi and signal Ei 


(3.5) 
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ZaP Plasma Emission Spectra 



Figure 3.1: ZaP line radiation. B-III (206.6 nm), C-III (229.7 nm), and N-III (409.7 
and 413.4 nm) represent the greatest contributors of line radiation in the ZaP plasma. 


divided by the respective area under the signal curve. 

^ _ Ya=i ^XjXiEj 


(3.6) 


The MATLAB code centroid.m used for centroid calculation is included in Ap¬ 
pendix A.4. These centroid data are important since they provide a comparison 
between plasma radiation and current. Comparing these data to the appropriate 
X- or I/axis component of the center of current data from the magnetic probes also 
provides valuable insight into the characteristics of the ZaP Z-pinch. 
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3.3 Centroid Error Analysis 

While the calculation of the centroid is somewhat straight-forward, there are many 
uncertainties inherent to experimental analysis. Standard deviation measurements 
indicate the degree of accuracy of the centroid calculation. Standard deviation (a) 
is a measurement of dispersion of the centroid due to error and is governed by the 
following equation [2], 


a = 


{EiXiAxi 



1 

2 


(3.7) 


Ei, Axi, and Xi are the measured signal, spacing, and position of each chord. 6Ei 
and 6xi represent the uncertainty of these signals and locations. In this case, the well 
known chord locations and spacings limit the relative error contribution from these 
terms to under 1%. The signal variation term however, contributes a higher error 
of ~ 5%. The error from this term is especially apparent for the flat profiles seen 
during formation and post-quiescent periods in the pinch, where variations between 
low Ei values are much more noticeable. Measurements of the centroid standard 
deviation determined by the MATLAB code centroiderror .m listed in Appendix A.5 
are plotted as error bars in Figure 3.2. 

Because centroid location plays an important role in future data analysis, it is im¬ 
portant to remove data points where there are noticeable differences between emission 
and the calculated centroid. The centroid. m algorithm removes data points from flat 
profiles and where the calculated centroid is beyond the outer channels of the array. 
If the maximum signal value of the array is less than 3 times the average of the two 
outside channels, the data are thrown out. This condition comes from comparing the 
movement of the emission profile to the calculated centroid, as shown in Figure 3.3. 
Further results for center and emissivity comparisons are listed in Section 5.3 and 5.4. 



1.5 



Time (|is) 

Figure 3.2: Centroid of ZaP emissivity and standard deviation. The distances from 
each centroid location indicated by the error bars represents the standard deviation 
given by Equation 3.7, giving upper and lower bounds to the emissivity location. The 
horizontal dashed lines represent the center of each channel along the PDA viewing 
axis. For this experiment, the standard deviation of the emissivity centroid was 
~ 1 channel spacing during quiescence (shown between 39.2 and 75.4 /rs indicated by 
the vertical blue dashed lines). When the emission prohle flattens out later in time, 
the standard deviation of the centroid increases to ~ 2 channel spacings. 


3.4 Chord-Integrated Emissivity vs. Density 

Chord-integrated emissivity from the photodiodes is also useful for comparison to the 
Z-pinch electron density. The radiation per unit volume for each unit of frequency is 
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Shot#50316017 



Time (jas) 

Figure 3.3: Truncation criteria for centroid data due to emission movement and flat 
profiles. The solid line represents the complete centroid data for this given emission. 
When the prohle at a given time step is not sufficiently peaked or is outside the 
viewing area of the array, the data are removed. Only the remaining points, shown 
in red, are used for comparison and further data analysis. The emission prohle shown 
here is saturated at 35% of the maximum value for viewing purposes. 


where P is power per electron and / (v) the electron distribution. If / (v) is isotropic, 
emissivity in terms of electron and ion density dehned by the term j [u) is [13] 


4:71 j (z/) = UeUiZ^ff 


3271 ^ 


2m 


1/2 




(3.9) 


^dvreo/ 

This equation represents the spectral power emitted over a full spherical solid angle 
per unit of frequency per unit volume. Ug and n* represent number densities, Zeg 


the effective ion charge, e the electron charge in Coulombs, eo the permittivity of free 
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space, c the speed of light, T the temperature at thermal equilibrium, hu the quantum 
energy, and g the Maxwell-averaged Gaunt factor. Taking the emissivity over areas of 
low line radiation and low photon energy reduces this equation to simply the number 
densities multiplied by a constant so that E (x n1 [12], Therefore, by viewing the 
wavelength region of the plasma emission with only bremsstrahlung, it is possible to 
determine a general density prohle of the pinch. The best spectra determined from 
Figure 3.1 for this condition fell between 510-580 nm. 

Unfortunately, these data are of little use for a conclusive calibration of the inter¬ 
ferometer using the PDA, as they view different volumes of the plasma and therefore 
cannot determine a correct scaling factor. However, by assuming axisymmetry of the 
plasma, it is possible to make general comparisons between the two measurements. 
Results for density calculations are listed in Section 5.2. 

3.5 Radial Profile Deconvolution 

One of the more exciting uses of chord-integrated data is for tomographic deconvolu¬ 
tion of linear chord data into a local radial prohle. By using one PDA with a linear 
array of photodiodes, it is possible to model the plasma as a combination of indepen¬ 
dent and piecewise constant cylindrical shells. This is accomplished using the well 
known Abel inversion technique [1]. 

3.6 Abel Inversion and One Dimensional Theory 

The emissivity viewed by each chord is merely the integration of all the radial prohles 
e (r) over the viewed length x, or E (y) = J e (r) dx. Taking advantage of the geometry 
of the path lengths through the plasma allows this integral to reduce to only chord y 
and radial r geometry. 



E{y) = 2 


(3.10) 
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Inverting this relationship to solve for e (r) requires determining the change in chord- 
integrated emissivity with respect to y, resulting in the following equation, 


e[r) = - 

71 


E'M 


Vv 


2 _ ^2 


dy. 


( 3 . 11 ) 


As there is only one viewing array along a single axis, it is important to note that an 
Abel prohle assumes that the plasma exhibits circular symmetry. 

With multiple viewing chords and discrete data, it is easier to solve the Abel 
inversion through the use of linear algebra and a numerical methods program such as 
MATLAB. The Abel inversion follows the general equation Le = E, where L is a path 
length matrix, e the radial emissivity, and E the chord-integrated emissivity. This 
method assumes that the radial prohle is piecewise constant over each shell. Figure 3.4 
shows a two chord example with Equation 3.12 as its matrix representation. 



(3.12) 


In MATLAB, the inversion takes advantage of the built-in least squares solver to 
increase computational speed. While the path length matrix in the one-dimensional 
Abel case is square, and therefore exactly determined, a least squares data-htting 
routine proved advantageous in the overdetermined two-dimensional deconvolution 
discussed later in this chapter. In the square matrix case however, this solver computes 
the solution using Gaussian elimination. A listing of the MATLAB routine abell.m 
is included in Appendix A.6. 


3.7 Simulated Data and Error Analysis for 1-D Radial Profiles 

Discrepancies are expected between any numerical solution and the exact analytic so¬ 
lution. In an attempt to reach a convergent solution to a linear problem, assumptions 
and approximations are always made [20]. Geometric considerations for the photodi¬ 
ode array, such as chord spacing, also play a factor in the accuracy of the results. In 



24 


y 



Figure 3.4: Abel inversion path length geometry 


order to get a better understanding of these influences on the Abel algorithm, it is 
important to check the deconvolved data against a known prohle. 

The MATLAB Abel inversion algorithm used to deconvolve the PDA chord- 
integrated data takes advantage of the centroid hnder used in Section 3.2 to divide 
the array prohle into left and right sections. Chord-integrated data are determined 
from known radial prohles, and a path length matrix is created from the known chord 
spacing. These data are then inverted and solved using Gaussian elimination. Peaked, 
hollow, centered, and offset prohles were all considered in this analysis. Table 3.7 is a 
list of several diherent synthetic prohles and their governing equations. These prohles 
correspond to the inversions shown in Figure 3.5. 

Increasing the number of chords viewing the plasma decreases the variation be- 
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(a) 



(c) 



(e) 



(g) 



Radius (cm) 


(b) 



(d) 



(f) 




Radius (cm) 


Figure 3.5: Synthetic Abel profiles for 16 and 32 chords. These correspond to the 
equations in Table 3.7. The black line represents the synthetic profile, the blue line 
the chord-integrated profile, and the red line the returned Abel profile. The points 
on these lines represent the chord locations. 
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Table 3.1: Synthetic profiles and governing equations for 16 and 32 chord examples 
used in Figure 3.5. Chords on each end of the prohle are set to zero in order to force 
the inverted prohle to zero at the edges. 


Prohle 

Equation 

Flat prohle 

(Fig. 3.7(a) and 3.5(b)) 

e (r) = 1 

r < 1 

Linear prohle 

(Fig. 3.5(c) and 3.5(d)) 

e (r) = 1 — r 

r < 1 

Parabolic prohle 

(Fig. 3.5(e) and 3.5(f)) 

e (r) = 1 — 

r < 1 



( ^.3 

for 0 < r < 0.5 

Hollow prohle 

(Fig. 3.5(g) and 3.5(h)) 

e{r) = < 





for 0.5 < r < 1.0 


tween the initial and deconvolved prohles. Reducing the separation between the 
centroid and chords immediately to its right and left also improved the continuity of 
the numerical routine between the two sides. 


3.8 Centering Errors of the Abel Inversion 

While the Abel prohles shown in Figure 3.5 are all properly centered on the given 
symmetric data, this luxury is unfortunately not always possible with a moving Z- 
pinch plasma. Although the inversion still assumes radial symmetry, any discrepancies 
between the center of the inversion shells and the center of the plasma prohle can have 
signihcant errors. For the abell.m inversion program, the inversion center is set to 
the centroid location determined by the program centroid.m. Figure 3.6 shows the 
discontinuity at the center of the Abel prohle caused by ohsetting the shell locations 
by 1 chord spacing. The prohle is given by the equation e (r) = 1 — r^, also shown 
correctly centered in Figure 3.5(f). 
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Figure 3.6: A 32 chord, offset synthetic profile given by the equation e (r) = 1 — 
The black line is the synthetic radial profile centered on the black dashed line. The 
blue line is the chord-integrated data taken from the radial profile. The red line 
indicates the radial Abel inversion profile with the shell centers set at 0 cm on the 
x-axis. 


3.9 Two Dimensional Tomographic Deconvolution 

While the Abel inversion technique is useful for resolving profiles with radial sym¬ 
metry, it cannot accurately model those with azimuthal variations. Inverting a non- 
symmetric profile can lead to discontinuities between the left and right sides. Rotating 
the viewing array into different axes also reveals changing radial profiles. Inverting 
chord-integrated data at cross-sections where the profile is not symmetric also leads 
to discontinuities between the left and right side sections, as shown in Figure 3.7. 

Simultaneously viewing the plasma in different orientations using multiple arrays 
allow for azimuthal variations to be considered in the deconvolution profile. This 
allows for better resolution of the emissivity of the plasma relating to the mode 
instabilities listed in Section 2.2. Fitting the chord-integrated data of two photodiode 
arrays to superimposed Fourier component profiles and solving a similar numerical 
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(a) Azimuthally varying hollow profile 




(b) Y-axis cross-section 


(c) X-axis cross-section 


Figure 3.7: (a) Abel inverted cross-sections of an azimuthally varying profile, (b) is 
the inversion taken along the x-axis, and (c) along the t/-axis, where the black line 
is the synthetic prohle, the blue line the chord-integrated data, and the red line the 
inverted prohle. The dots represent the discrete data points. Note the discontinuity in 
(c) resulting from the Abel algorithm attempting to follow the non-symmetric chordal 
data across the left and right sections. The hollow prohle in this example follows the 

, , , , , , ^ , , . [ for 0 < r < 0.5 

equation e (r) = en (r) -|- O.Sen (r) cos6', where en (r) = < o 

^ ^ ^ ^ ^ ^ ^ ^ 1 (1 - r)^ for 0.5 < r < 1.0 
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linear algebra routine as the one-dimensional case allows for the generation of a shell 
prohle of the Z-pinch varying both radially and azimuthally. 

3.10 Fourier Component Theory 

Solving a two dimensional plasma prohle involves using a Fourier series expansion 
method developed by Cormack [4, 5]. A radially and azimuthally varying plasma 
prohle is represented in polar coordinates by g (r, 6). For a chord viewing a circular 
plasma through a length L, the line integral of the prohle is simply 

= [ g{r,e)ds. (3.13) 

p represents the perpendicular distance from the geometric origin of the prohle to the 
line L, and 0 is the angle from the positive x-axis to the line measuring p. Because 
of the viewing geometry for the ZaP experiment and the restricted use of only two 
PDAs, the value of p is equal to each array’s respective chord location along the 
viewing axis. 0 is also limited to a value of 0° or 90°, again depending on the viewing 
axis. In addition, if g is assumed to be piecewise constant, this expands in Fourier 
series to 

1 

( 3 . 14 ) 

n=0 

where gn is the piecewise constant prohle for shell m. Without this is simply 
the Abel prohle described in Section 3.6. Expanding the complex into sine and 
cosine components dehned by e*"® = cos nd + i sin nd results in 

g (r, 6*) = go (r^) + gi (r^) [cos6' ^sin^]. (3.15) 

This equation represents a non-symmetric radial shell prohle by combining the de¬ 
convolved Abel prohle with one multiplied by cosine and sine functions. 

As with the Abel inversion algorithm, the two-dimensional prohle lends itself to 
matrix computation and linear algebra. Rewriting equation 3.15 so that each shell 



30 


profile e (r) = cq + ei cos 6 ' + 62 sin 6 *, the relationship between chord-integrated data 
and radial data takes the form of E = Me, where E are chord-integrated values, 
M are the combined path length and sine and cosine components, and e are the cq. 
Cl, and 62 values for each shell. A criterion that at least three chords pass through 
the outer shell is used to determine the maximum number of shells deconvolved for 
a given prohle. Smaller shells are added inside of this one at a distance equal to the 
average chord spacing. This ensures that the number of chords provided by the PDAs 
will always be greater than or equal to the unknown cq, ci, and 62 prohles for all the 
shells, or rather that it is either exactly or overdetermined. An elementary one shell 
example is given in Figure 3.8, along with the corresponding matrix represenation in 
Equation 3.16, where 6 is the angle from the positive x-axis to the midpoint of each 
pathlength element. 


Ell 


Ei2 


E 21 


E 22 



2Lii 

2Ti2 

2 T 21 

2T22 


Lii cos 011 -|- Lii cos (—0ii) 
Li2 cos 012 -I- Li2 cos (-012) 
L 21 cos (-022) + T2I COS (-02l) 
L 22 COS 022 + L 22 COS 021 


Lii sin 011 -I- Lii sin (-0ii) 
Li 2 sin 012 -I- Li 2 sin (-0i2) 
L21 sin (-022) + L 21 sin (-02i) 
L 22 sin 022 -I- L 22 sin 02i 



(3.16) 


Unlike the 1-D Abel case, this method solves for all the chords from both the top 
and side viewing arrays in one single path-length matrix instead of dividing them by 
their respective centroid locations. Because the path-length matrix M can be over¬ 
determined and therefore rectangular, the solving algorithm requires a least squares 
fit routine to calculate the best inverted prohle. The MATLAB routines TDerror.m 
and bestf it .m make use of QR decomposition with pivoting to converge to the best 
solution. 9 values in this matrix represent the angle from the positive x-axis of the hrst 
PDA to the midpoint of the respective chord through an individual shell. This small 
angle approximation becomes less pronounced with an increasing number of chords. A 
listing of these two MATLAB routines and the overall Fourier deconvolution algorithm 
FULL2d.m are given in Appendix B. 
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Figure 3.8: 2-D Fourier deconvolution path length and angle geometry. Line segment 
midpoints are used to determine 6 angles for cosine and sine components. 


3.11 Simulated Data and Error Analysis for 2-D Radial Profiles 

Comparisons of the deconvolved data to synthetic radial profiles are just as important 
in understanding the limitations of the 2-D Fourier deconvolution method as they were 
for the 1-D Abel inversion algorithm. Validating the code involved using data from 
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azimuthally varying profiles, both centered and offset with respect to the photodiode 
arrays. Results exposed several different factors that play an important role in the 
amount of error generated in the 2-D numerical method. 

In order to minimize the approximations and converge to a solution, the least 
squares £t solver attempts to find the minimum norm, ||a;|| 2 , with no rank deficiency. 
A rank deficient matrix indicates that several channels are not linearly independent, 
meaning that there are not enough data points to resolve unknowns in some of the 
shell profile. This occurs primarily from not having enough independent chords for 
either the outer or inner shell profiles. Minimizing these conditions involved setting 
the maximum shell radius to the criteria set in Section 3.10 and by guaranteeing a 
minimum of four chords traverse the inner shell. Although only three chords are 
required to solve for the three unknown eo,i ,2 values, when a chord passes very near 
the origin of the inner shell, the larger jump in path-length sizes and the proximity of 
their 6 angles to 0 or 90° can skew the data. Increasing the inner shell radius to allow 
an extra chord to pass through reduces the chances that this condition will exist. 

Correct geometric centering of the inverted profile on the plasma center is also 
a serious concern in error analysis. Since path-lengths passing through the origin 
of the radial model can distort the profile of the inner shell, the algorithm forces a 
minimum perpendicular separation between any given chord and the chosen center. 
Unfortunately, since the entire shell model geometry revolves around this point, any 
offset in the routine center can also have disastrous consequences and result in a 
deconvolved profile very different from the synthetic contour. A minimum error finder 
determines the best profile center by comparing the given and the calculated chord- 
integrated data and finding the minimum location, is a measurement of the 
goodness of fit of least squares data, and is governed by the following equation [2], 

= ^[^[yi-y{xi)f^ . (3.17) 

The standard deviation (a) squared term is the variance, or expected error in the 
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Table 3.2: Synthetic profiles for 2-D Fourier analysis used in Figures 3.9(f)~3.11(f). 
An azimuthal variation is added to e (r) by sine and cosine components 30% the size 
of the axisymmetric data. 


Prohle 

e (r) 

= e -f 0.3e cos 6 * -|- 0.3e sin 6 

Parabolic 

(Fig. 3.9(a), 3.10(a), 3.11(a)) 

e = 1 

_ ^2 

r < 1 




^3 

for 0 < r < 0.5 

Hollow 

(Fig. 3.9(d), 3.10(d), 3.11(d)) 

e = < 






(l-rf 

\ 

for 0.5 < r < 1.0 


signal of each channel, yi represents the measured signal and y (xi) is the returned 
calculated signal. When the htting routine is a good approximation of the given syn¬ 
thetic data, the value of is on the order of 1. Contour plots of several azimuthally 
varying prohles and their respective data are shown in Figure 3.9. Table 3.11 is 
a list of the equations used to generate the synthetic prohles. These data show that 
the width of the contour of = 1 is always much less than the spacing between 
the viewing chords, meaning that the hnal position of the prohle center chosen by 
the algorithm is very well known compared to the chord spacings. The MATLAB 
algorithm FULL2d.m and subroutines TDerror .m and bestf it .m calculate the decon¬ 
volved 2-D prohle at the minimum position. Keeping track of location over time 
also allows the algorithm to track the plasma prohle movement oh the ; 2 -axis, as shown 
in Figure 3.10. 

The TDerror .m and bestf it .m routines also compensate for cases where the large 
variation in sine and cosine contributions lead to incorrect prohles. This involves 
constraining the value of the axisymmetric cq term by the ei and 62 data according 
to the following equation, 

|eo| ^ + ei. (3.18) 

Without this constraint, the sine and cosine components of the deconvolution can 
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cause the routine to return a negative emissivity profile. These cq, ei, and 62 values 
are used to determine the calculated e (r) values. Multiplying e (r) and the path- 
length matrix M produces the calculated signals y {xi) used in minimization. 

In addition, random noise is added to the algorithm to test the effects of digitizer 
noise in prohle reconstruction. This is accomplished by adding random noise with a 
known magnitude to the chord-integrated data. The a* term in Equation 3.17 is set 
to this magnitude in order to match the known error. During experimentation, Uj is 
determined by calculating the bit noise of each digitizer channel. Figure 3.11 shows 
the effects of adding a 5% error to the input signal. 
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Figure 3.10: 2-D synthetic and deconvolved profiles with x and y = -0.5 cm offset. The algorithm for the circular 
prohle set the reconstruction center at x and y = -0.43 cm. The deconvolved hollow prohle is centered at x and y = 
-0.53 cm. Both prohles are still within a single chord spacing of the synthetic data. 
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Figure 3.11: 2-D synthetic and deconvolved profiles with added noise equivalent to 5% of the maximum signal. The 
noise has an average of zero and a standard deviation of 1. Again, the goodness of £t determined by = 1 varies by 
less than the chord spacing. 
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Chapter 4 

PHOTODIODE ARRAY DESIGN AND CALIBRATION 

Using a photodiode to detect plasma emissivity is a straightforward process that 
simply involves pointing it at the light source and measuring the resulting current. 
Unfortunately, taking useful tomographic data requires a great deal more design con¬ 
sideration and calibration. An array of photodiodes collects light from different re¬ 
gions of the plasma. Collimation of these photodiodes reduces the visible solid angle 
and provides meaningful data based upon the array conhguration. This requires a 
well known array geometry and dehned viewing angles that remove reflected light. 
The vacuum tank window design and inner electrode openings further restrict the 
area of the measured plasma. Components of the overall photodiode array used in 
experimentation are given in Figure 4.1. Knowledge of the spectral response of the 
photodiode chip, individual signal amplihcation and attenuation, frequency response 
and noise hltering are all required to accurately interpret the resulting data. 

4.1 Photodiode Chip Design, Spectral Response, and Biasing Board 

As previously mentioned in Section 2.1, an aligned array of photodiodes allows for 
the collection and generation of an emissivity prohle of the Z-pinch plasma. The 
ZaP experiment makes use of four, 16-element silicon photodiode arrays made by the 
Hamamatsu Corporation. The S5668-02 series proved ideal in providing a quick signal 
rise time and linear array configuration for tomographic analysis over the visible light 
spectrum. Each PDA uses two aligned chips to create a 32 channel linear array. Since 
the photodiodes all have a common cathode, they produce a small negative closed loop 
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TANK CONNECTION gS>4.62 .n diameter 

FLANGE 0.20 in thickness 

0.332 in bolt holes 
2.50 X 0.75 in opening 


Figure 4.1: Components of the photodiode array 
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current when exposed to a light source. Each active element has a 1.175 mm width 
and 2 mm height along the viewing axis. The mounting PC board is 25 mm (1 in) 
in length, allowing for easy mounting and alignment. A schematic of the photodiode 
chip geometry is shown in Figure 4.2 [10]. Table 4.1 lists typical electrical and optical 
characteristics of the photodiode chip. 


NDTEi DIMS IN 

Conrion Cathode Pin 



1,0 


FeNiCo PIN 


-Glass Epoxy Board 
G-10 (Black) 

Epoxy Resin 


3,5 


Figure 4.2: Hamamatsu S5668-02 silicon photodiode chip. Dimensions are in milli¬ 
meters. 


As shown in Table 4.1, the S5668-02 series photodiodes have wavelength detection 
limits between 320-1100 nm, with a peak sensitivity of 960 nm. The photocurrent 
produced from incident light between these boundaries varies according to wavelength 
along a spectral response curve, shown in Figure 4.3 [10]. The higher order lines seen 
in the spectra plots represent the light scattered in frequency space by the diffraction 
grating of the ICCD spectrometer. They do not contribute to the actual amount 
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Table 4.1: Important electrical and optical characteristics for the Hamamatsu S5668- 
02 silicon photodiode chip. 


Parameter 

Symbol 


Unit 

Spectral Response Range 

A 

320-1100 

nm 

Peak Sensitivity Wavelength 

Xp 

960 

nm 

Photo Sensitivity 

s 

See Figure 4.3(a) for curve 

A/W 

Rise Time 

tr 

0.1 

fiS 


of line radiation in the plasma and can therefore be ignored. As previously seen 
in Figure 3.1, B-III, C-III, and N-III ions produce a great deal of line radiation. 
Fitting the response curve to the plasma radiation plot reveals that B-III and C- 
III are completely hltered out by the photodiode chips. N-III remains a signihcant 
contributor to emissivity for these pulses, however higher wavelength line radiation 
now make up a larger proportion of the light detected. 

Since plasma instabilities in the ZaP sheared flow Z-pinch occur on the order of 
microseconds or faster, it is imperative to increase the response time of the photodiode 
circuit to this same order of magnitude. A biasing circuit running off of a 9 V battery 
accomplishes this by introducing a reverse bias current on the photodiodes. This 
current replenishes the electron hole pairs in the silicon crystal depletion layer between 
the p and n sections of the photodiode, thereby reducing the signal current delay. This 
board mounts directly onto the back of the photodiode chips, as shown in Figure 4.4. 


4.2 Collimator Design 

Alignment and collimation of the photodiode chips is accomplished using a DuPont 
plastic polymer called Delrin. Delrin proved advantageous because of its durability. 
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ZaP Plasma Emission Spectra 



Wavelength X (nm) 

(a) 


ZaP Plasma Emission Spectra 
as seen by the Photodiodes 



Figure 4.3: Effects of photodiode sensitivity on viewed emission spectra, (a) is the full 
emission spectra with the the photodiode crystal photo sensitivity curve overplotted, 
(b) is the ZaP spectra as seen by the photodiodes due to the response curve. N-III 
remains the top supplier of line radiation. Higher order lines seen in the spectrometry 
measurement may be ignored. 
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Figure 4.4: PDA biasing board mounted on the photodiode chips. Each set of twisted 
pair wire connects to one of the thirty-two channels from the photodiodes. 


ease to machine, smooth surfaces, and its availability in a black opaque color that 
would minimize reflected light. Four pieces of Delrin plastic were used in the con¬ 
struction of each of the two collimators, and most of the geometric restrictions were 
a direct result of limits from the window mounting on the ZaP vacuum tank and 2 in 
diameter outer electrode view ports. A previous sixteen element photodiode array 
device used on the ZaP experiment offered a template for the base design of the cur¬ 
rent arrays. An exploded view of the current collimator design is given in Figure 4.5 
and a color enhanced schematic of the combined assembly is shown in Figure 4.6. 

The photodiode chips are housed in a recess in a 3.6 x 1.75 x 0.25 in Delrin cradle 
piece. The sides of the Delrin opening provide horizontal alignment of the chips 
above thirty-two 0.032 in diameter holes. Each hole is centered 0.062 in apart along 
the horizontal axis with tolerance levels minimized to approximately 0.001 in through 
the use of a digital display endmill. The chips are secured by backing socket pins and 
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Figure 4.5: Unassembled photodiode array collimator 
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Figure 4.6: Assembled and aligned photodiode array collimator 
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the biasing board, which attaches to the Delrin cradle using four 8-32 screws. The 
cradle assembly sits atop two interconnecting side pieces, a large smooth side piece 
and the slot piece. The slot piece has thirty-two 0.040 in wide slots with 0.022 in 
ridges and is 3.6 x 0.91 x 2.78 in in size. Alignment between the cradle holes and the 
slot piece occurs through the use of two 0.187 in alignment pins, secured in specially 
drilled holes. Four 8-32 bolts connect the sides together like puzzle pieces. A small 
0.020 in thick opaque rubber piece sits between the smooth side and the slots. The 
compression force created by the bolts ensures that this rubber piece will close any 
gaps between neighboring slots, reducing the chance that light will “leak” from one 
channel to another. The side pieces combine to produce the same length and width 
prohle as the cradle piece. A “top plate” with the same dimensions, alignment pins, 
and hole conhguration as the cradle attaches to the other end of the side pieces 
to create a rectangular box collimator. A smaller recess was cut into this piece to 
allow for easier drilling of the collimation holes, and to minimize difference in pitch 
caused by the drill press. In all, the four pieces combine to offer approximately 2.9 in 
of collimator length from the top plate opening to the photodiode crystal. Overall 
design tolerance for the assembly was approximately 0.003 in. 


Since the signal detected by the PDA is directly proportional to the amount of 
light striking each photodiode crystal, it is important to know the viewing solid angle 
of each channel. Although the collimator design highly parallelizes the light, by the 
time the light cones reach the center of the tank (7 in from the vacuum chamber 
window), the viewing area created by the 0.032 in (0.08 cm) diameter collimator hole 
balloons to a diameter of 0.19 in (0.49 cm). As shown in Figure 4.7, this creates a 
total viewing area of 2.11 x 0.19 in (5.37 x 0.49 cm) where adjacent channels overlap. 
This overlap acts to smooth the chord-integrated profile recorded by the PDA. Using 
fewer channels on the PDA allows for the generation of a more discrete prohle. 
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Figure 4.7: Collimator viewing areas of the ZaP plasma. All dimensions are given in 
centimeters, (a) shows the collimator hole openings and (b-d) the channel viewing 
angles of a 1 cm radius plasma 7 in from the collimator circular connection flange, 
(b) shows 8 equally spaced channels, (c) 16 channels, and (d) all 32 channels. 


4.3 Amplifier Board Design 

While the collimator excels at highly parallelizing the viewed light, it has the unfor¬ 
tunate consequence of signihcantly reducing the photodiode output current. Back¬ 
ground noise is therefore much more signihcant, and special care must be taken to 
shield collection circuits from the lab environment and hlter out unwanted noise. The 
amplihers also match the output impedance of the photodiodes to the 50 input of 
the digitizers. 

Each photodiode in the array has an amplihcation circuit consisting of two op- 
amps; a transimpedance ampliher and a signal buffer. Figure 4.8 shows the electrical 
schematic for a single ampliher. The AD797 op-amp provides signal amplihcation 
with very little noise distortion. Unfortunately, it also has a small 50 mA peak output 
current. A subsequent LM7171 buher increases this output current to a usable level 
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Trimmer Capacitor 


1.5-lOp 



Figure 4.8: Photodiode amplifier circuit diagram 


and provides voltage amplification of the final signal. The 47 hi resistor following the 
LM7171 prevents transmission reflections from the cables in the lab and protects the 
amplifier in the event of an electrical short. The only variable component of the design 
is a 1.5-10 pF trimmer capacitor, which allows for frequency filtering by correcting 
underdamped or overdamped signal outputs. 

High background noise levels in the lab caused by the main bank power supply 
and plasma current interruption during post-quiescence are of great concern to the 
photodiode array design as they can interfere with the collection of any meaningful 
data. The amplifiers filter out this noise by doubling as low pass frequency filters. 
Frequency response tests on the amplifier, as shown in Figure 4.9, indicate that they 
successfully remove high frequency noise from data with a frequency below 2 MHz 
without negative gain. Because the frequencies of interest for the ZaP experiment are 
well below this 2 MHz corner, any remaining high frequency noise can be neglected by 
smoothing the data. The amplifiers, collimator cradle, biasing board, and photodiode 
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PDA Frequency Response 



Figure 4.9: Photodiode amplifier frequency response adjusted to DC gain = 1. The 
black points represent signal gain at given input frequencies. The red line is a rough 
trendline of these points. These data show that the amplifiers can filter out high 
frequency noise without signal degradation below an input frequency of 2 MHz. 


chips are all shielded from the plasma magnetic field by a thin aluminum box, as shown 
in Figure 4.10. The output signals for the thirty-two channels are divided between 
three DB-25 connectors on a removable side panel of the box. In order to fit everything 
inside the box, the amplifier circuits are mounted onto PC boards and stacked using 
8-32 threaded rod and insulating spacers. The sixteen odd channels and sixteen even 
channels are divided into two groups of four boards, each receiving a separate ±12 V 
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Collimator, Biasing board, mid Photodiode chips 


Figure 4.10: Stacked amplifier boards inside the shielding aluminum box 


power supply required by the op-amps. This division of the input power is important 
because of the high supply current required by the ampliher boards. Each side of the 
AD797 and LM7171 chips have peak input currents of 10.5 and 9.5 mA respectively. 
The thirty-two channels add up to an overall peak supply current requirement of 
1.28 A. Dividing the input current into two equal parallel supplies reduces the risk 
of overheating any of the circuit components. Thermocouple tests conducted on the 
PDA’s during data collection showed that element temperature levels remained well 
below the 85° C (185° F) maximum operating temperature using this power supply 
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configuration. The external power is also divided between the DB-25 connectors 
and twisted pair cable running out of the laboratory. A picture of the overall PDA 
assembly mounted on the ZaP Experiement is shown in Figure 4.11. 


4.4 Photodiode Array Channel Calibration 

Due to variations in the photodiode chip, collimator, biasing and ampliher boards, 
calibration factors are required to correlate the channels on both arrays to one relative 
standard. This is accomplished by relating the maximum signals detected by each 
channel from a constant light source. During the testing and calibration phase, the 
photodiode array was secured on a translational optics table such that the moving axis 
was aligned with the photodiode viewing axis. A Strobotac® 1531-AB series strobe 
light made by the General Radio Company with a flash frequency of approximately 
167 Hz was placed behind a frosted glass window and an opaque cardboard covering 
with a narrow 1.0 in slit window. The frosted glass and slit window ensured that only 
a small band of uniform light was visible to the PDA. The strobe light was placed 
on the translational table approximately 7 in from the collimator opening, roughly 
the same distance it would be from the center of the plasma when in use on the 
experiment. Figure 4.12 shows a picture of the channel calibration setup. An external 
±12 V supply was applied to power the PDA ampliher boards. An oscilloscope was 
attached to the DB-25 connectors to measure channel output signals as the strobe 
light translates along the optics table. The peak to peak output signal is recorded for 
each channel when the slit window was aligned directly over that particular collimator 
hole. The movement of the table also provides a rough measurement of the distances 
between the viewing angle centers of adjacent channels. The average distance of 
0.15-0.16 cm measured corresponds to the distances between ray traces in Figure 4.7. 
Figure 4.13 shows an example of the measured strobe light signal and translation 
as seen by the oscilloscope. This testing also provided a means to check ampliher 
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Figure 4.11: Photodiode array mounted on the ZaP experiment 
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Figure 4.12: Channel calibration setup 


(a) 



(b) 



Figure 4.13: Channel calibration measurements as seen by the oscilliscope. Four 
output channels are displayed in each plot. In Figure (a), the strobe light slit is 
aligned over the collimator hole for ch. 1. In Figure (b), the strobe light is translated 
so that the collimator hole for ch. 4 is now aligned with the slit. The peak-to-peak 
voltage measured for each channel is used to create calibration factors. 
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damping levels and circuit connection problems inside the amplifier box. Over and 
underdamped signals were adjusted to corrected levels using the trimmer capacitors 
on the amplifier boards, and any noise spikes or loose connections were investigated 
and fixed before recalibrating that channel. Figure 4.14 illustrates several different 
problems encountered during channel calibration. 



(b) 




Figure 4.14: Amplifier board errors. Figures (a) and (b) show a fiatline and internal 
noise signal respectively. These are indicative of a loose solder connection on that 
particular channel. Figure (c) shows an underdamped amplifier response. Adjust¬ 
ing the trimmer capacitor shown in Figure 4.10 fixes this problem, however a new 
calibration factor must then be taken for that channel. 
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The biasing and amplifier boards were also tested independent of the collimator 
for cross talk and frequency response characteristics. With the socket removed from 
the photodiode chips, a square signal from a waveform generator was supplied to 
the circuit as shown in Figure 4.15. The frequency response of the amplifier circuits 
calculated from this test is approximately 8 MHz, well above the frequency of interest 
for this experiment. A second oscilloscope channel measured the effects of cross talk 
between adjacent channels on the amplifier boards. The noise measured by the second 
channel was under 1% of the input square waveform and was therefore negligible. 



Figure 4.15: Amplifier board frequency response and channel cross talk measurements. 
Channel 1 measures a rise time of 127.4 ns associated with the square waveform input. 
This corresponds to a frequency of 7.85 MHz. Channel 2 has no direct input signal, 
but it measures a small peak to peak voltage induced from the adjacent channel’s 
waveform. 
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4.5 Data Collection Components and Calibration 

For experimental operation, the PDA collimator and amplifier box are mounted on 
the ZaP vacuum tank’s 4-5/8 in confiats with 3 in diameter fused silica view-ports 
using a circular aluminum flange with a wide slit opening in the center, as shown 
in Figure 4.16. Section 4.6 includes a more detailed discussion of the tank windows 


PDA 2 



Z = 0 Plane 


Faraday Cage 


+ Y 


Figure 4.16: Photodiode arrays mounted on the ZaP experiment. Note that the 
vacuum tank extension is not shown on the schematic. All experimentation was 
conducted at the z = 0 plane. 


and inner electrode openings. From the DB-25 connectors on the amplifier box, 
each output signal travels to analog digitizers along 50 ft of 7 strand, 20 gauge, 32 
twisted pair wire with both individual and overall shielding and grounding. To avoid 
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generating gronnd loops, the cable is gronnded on one side to the DB-25 connectors 
on the PDA box and again at the entrance of the Faraday cage, or screen room 
containing the digitizers. The signal wire is converted from the twisted pair wire to 
Ethernet cable nsing RJ-45 jacks. These are then converted to BNC pings monnted 
on an alnminnm box at the top of the screen room, as shown in Fignre 4.17. BNC 
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Fignre 4.17: Screen room box connections. The box contains the gronnding wire to 
the screen room, conversion from the 32 twisted pair cable to RJ-45 jacks and ethernet 
to BNC plugs. 





58 


to LEMO connector cables connect this box to two, 16 channel Joerger Model TR 
analog digitizers. Each 12 bit digitizer operates at 40 MHz and has a ±2 V input 
range. Unfortunately, the voltage magnification by the amplifier boards requires the 
use of passive signal attenuation to reduce the digitizer input signal to the range 
indicated. 10 and 14 dB attenuators placed in-line between the BNC pins on the 
screen room box and the BNC to LEMO cables as shown in Figure 4.17 prevent the 
input signals from saturating the digitizers. 

Attenuation factors are generated by measuring a known waveform with a 1- 
2 MHz frequency through each attenuator using the circuit diagram in Figure 4.18. 
The attenuation value is simply the voltage out of the attenuator divided by the 



Figure 4.18: Attenuator calibration circuit schematic. In hgure (a), a known voltage 
Vo is applied in parallel through an attenuator and a 1 MVL resistor. Figure (b) is an 
enlarged view of the attenuator. 


voltage in, given by the equation 

Ui 


5OR3 


Vq 2i?i/?3 -|- 4- 50 (i?i -|- Rs) 

where resistor Ri = R 2 and follow the relationship 


-1 


/I 1 

50 = Ri+ { -^- 

\R3 Ri + 50 

These values are included in the calibration factors used for data analysis. 
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4.6 Viewing Constraints and Window Calibration 

Variations in the vacuum tank window and inner electrode opening size also affect the 
signal detected by the photodiodes. As shown in Figure 4.19, the two PDAs view the 
plasma through a 3.0 in diameter fused silica window and a 1.86 in diameter outer 
electrode hole. This small opening restricts the viewing axis of the array such that 
the 1st and 31st collimator holes clip the sides of the outer electrode, and channel 32 
does not directly view the opening. Therefore, the data taken by the arrays represent 
the complete observable area through the view-ports. 

The fused silica vacuum tank window also proved to be a significant contributor 
to variations in the collected data. During the course of experimentation, the plasma 
deposits tungsten metal and other various elements onto the view-ports. The accrue¬ 
ment of this material distorts the signal seen by the PDAs by diffusing the plasma 
emission through the window. Attempts to calibrate these distortions proved incon¬ 
clusive and led to the removal and sonic cleaning of the fused silica. The Strobotac® 
strobe light used for channel calibration was mounted adjacent to the view-port as 
shown in Figure 4.20 and set to high luminosity with a hash frequency of 5 Hz. A 
frosted glass window was also placed between the strobe and the view-port opening 
to achieve a more uniform intensity, and the Joerger digitizers collected data from a 
single light pulse by using an optical cable to trigger a measurement. Averaged plots 
of these data through the old and new window, as well as through no window at all 
are presented in Figure 4.21. 

The plots in this figure show significant diffusion of the strobe light through the 
old windows that would make any sort of relative comparison between the two PDA’s 
questionable. While the new window and no window data are similar from channel to 
channel, remaining variations in the channel calibration factors and frosted glass, non¬ 
uniformity of the source light, and misalignment of the strobe light and viewing PDA 
preclude this relationship from being used to generate window calibration factors. 



Figure 4.19: Viewing area of the Z-pinch seen by two PDAs at the z = 0 windows. 
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Frosted Glass Strobe Light 



Vacuum 
tank 
4-5/8 in 
conflat 
with 3 in 
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Stand 

Figure 4.20: Window calibration of the PDA on the ZaP tank. The photodiode array 
being tested is mounted on the opposite window viewing the strobe light. The fused 
silica window adjacent to the strobe is removed for the test, and data are taken for 
the dirty window, new window, and no window on the PDA side. 


Further attempts to use an opal glass lens to better diffuse the strobe light also 
proved inconclusive. In the end, shots taken with the clean windows are assumed to 
have no attenuation across all of the channels. 
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(a) PDA 1 Window Calibration 



-^PDA1 Old Window ♦PDAI No Window -^PDA1 New Window 


(b) PDA 2 Window Calibration 



Figure 4.21: Averaged frosted glass window calibration data for both PDAs. The 
frosted glass placed between the strobe and the tank as shown in Figure 4.20 produced 
a centered prohle instead of the desired nniform prohle. Note that the debris bnild- 
np on the old windows has signihcantly rednced and diffused the signal. The spike 
at channel 22 was the result of an improperly terminated digitizer channel. Other 
channel variations are primarily dne to ampliher calibration problems and variations 
in the frosted glass. Withont a trne nniform prohle, it is not possible to prodnce 
accnrate window calibration factors. 
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Chapter 5 

EXPERIMENTAL RESULTS 

After cleaning the vacuum tank windows, all emissivity data from the PDAs were 
collected at the z = 0 plane indicated in Figure 4.16. A bank charge of 5-9 kV was 
utilized for all shown data. The majority of the data listed in this chapter are from 
the 5 kV shots, where more useful data was available due to a longer quiescent period. 
Ha and bremsstrahlung hlters were also used to view different regions of the spectra. 
The bremsstrahlung hlter in particular allows for the comparison of PDA emissivity 
against chord-integrated electron density measured by the interferometer, as previ¬ 
ously discussed in Section 3.4. Comparisons between measured Imacon and PDA 
emissivity, and magnetic center of current and emission centroid are also determined 
from these chord-integrated data. 

An axisymmetric, piecewise constant radial prohle of the ZaP Flow Z-Pinch was 
determined from the 1-D Abel inversion technique described in Section 3.6. The 2-D 
Fourier deconvolution technique from Section 3.9 provides a more detailed view of the 
plasma, with both a radial and azimuthal piecewise constant prohle. This technique 
also tracks the movement of the prohle oh the 2 :-axis and can be used to investigate 
variations between the emissions of unhltered, Ha, and bremsstrahlung hltered light. 
Discussions of these results are included in this chapter. 

5.1 Chord-Integrated Results 

PDA emissivity data stored on the Joerger Model TR digitizers were run through the 
MATLAB codes using the theory described in Chapter 3. Each PDA used 16 channels 
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to collect data through the x- and y-ax.es of the tank at the z = 0 plane. The ma¬ 
jority of the channels are evenly spaced; however, problems resulted when variations 
between channel calibration factors resulted in signals with disproportionably larger 
or smaller gain than their neighbors. Studying the raw output data revealed where 
these problems were pronounced, and care was taken to take emission measurements 
with channels where these signal variations were minimal. Unfortunately, a bad at¬ 
tenuator on the screen-room box shown in Figure 4.17 resulted in the removal of one 
channel from the x-axis viewing PDA. 

Magnetic data taken for a 5 kV shot shows a dehned current sheet during plasma 
formation and an approximately 1 cm radius peaked emission prohle during the quies¬ 
cent period described by the magnetic probes. These results are shown in Figure 5.1. 
The unhltered emission also shows movement of the pinch along each independent 
viewing axis, as well as pronounced amplitude “flashes” during quiescence. 

Filtered data also provide useful information about the ZaP Flow Z-Pinch. Filter¬ 
ing the viewed spectra involved simple Kodak-Wratten 7^12 and 7^58 hlters. These 
him hlters are placed between the collimator top plate and the vacuum tank window. 
The line radiation detected by the photodiodes previously shown in Figure 4.3(b) 
shows that N-III is the primary contributor to emissivity. These lines are further 
reduced by the hlters according to the response curves shown in Figure 5.2. The 7)^12 
hlter, (a) and (b), only allows transmission of radiation above 510 nm, ehectively 
eliminating all major lines other than Ha- Ha is of interest because it represents the 
emission given oh due to collisions between neutral gas in the plasma with hotter ions. 
This sort of interaction occurs at the boundary of the pinch, therefore a emission hl- 
tered for Ha should have a hollow prohle at the pinch radius. The 4^58 is a band-pass 
hlter that only passes radiation between 500-580 nm. When combined with the 7)^12 
hlter, only radiation between 510-580 nm is detected by the photodiodes. This region 
of the ZaP plasma has very little line radiation, therefore combining the hlters allow 
the PDAs to measure primarily continuum bremsstrahlung radiation. 
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Normalized m = 1 Magnetic Mode Dispiacement 



Figure 5.1: Unfiltered chord-integrated results for 2 PDAs viewing a 5 kV bank voltage 
plasma. Peaked profiles with an approximately 1 cm radius are observed during the 
quiescent period described by the normalized m = 1 mode data and shown between 
the vertical dashed lines on all plots. Pinch movement along the viewing axis is also 
observed during this time, as indicated by the centroid position shown in red. During 
the formation period, the PDAs record the current sheet traversing the z = 0 plane, 
shown in this pinch at approximately 24-27 fis. 
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ZaP Plasma Emission Spectra 
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Figure 5.2: Filter curves and viewed spectra of the ZaP Flow Z-Pinch. (a) The blue 
line is the spectra seen by the photodiodes without hltering. The red line repre¬ 
sents the spectral response curve of the Kodak-Wratten ^12 hlter. (b) This hlter 
leaves Ha as the primary contributor to line radiation. (c,d) The response curve and 
viewed spectra for the ^^58 hlter still leaves some line radiation around 500 nm. (e) 
Combining the two hlters produces a narrower response curve between 510-580 nm, 
thereby removing almost all line radiation from the ZaP plasma, (f) This leaves only 
bremsstrahlung radiation as the primary emission seen by the PDAs. 
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Comparisons between 5 kV bank voltage shots for these hlters are shown in Fig¬ 
ure 5.3. Although the hlters signihcantly reduce the brightness of the measured light, 
especially in the bremsstrahlung case, the PDA is still able to detect a peaked prohle 
during quiescence. These data are therefore useful for further analysis later in this 
chapter. 

5.2 Electron Density Results and Comparisons 

Evaluation of the emission prohle from bremsstrahlung radiation and data from the 
interferometer shows that the electron density of the plasma is proportional to the 
square root of the emissivity. Interferometry data taken through the y-axis at z = 0 
prior to mounting the side-mounted PDA provides line-integrated electron density 
of the plasma along a chord through the center of the window at z = 0 and along 
a second chord 1 cm oh axis. These data show a peaked density plasma prohle 
during quiescence, with higher density measured along the center chord. Although 
the interferometer is removed to allow both PDAs to view the plasma at z = 0, 
relating 5 kV bank voltage shots with similar m = 1 magnetic mode data allows 
for general comparison between the interferometry and PDA signals. As described 
in Section 3.4, bremsstrahlung emissivity E oc n^. Figure 5.4 shows that the signal 
detected by PDA channels at similar centered and ohset locations with respect to the 
X- and y-axes also display this same peaked prohle. Since the PDA and interferometer 
are looking at diherent plasmas, amplitude variations in the two chords indicative of 
pinch movement are not similar; however, the signal at the center chord location 
generally remains higher than the ohset channel during quiescence. 

5.3 Visual Emission Results and Comparisons 

The Imacon streak photos previously discussed in Section 2.4 provide an independent 
verihcation of the photodiode emissivity. Since the time scale for the Imacon photo is 
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Figure 5.3: Chord-integrated profiles of (a) unfiltered, (b) Ha, and (c) bremsstrahlung 
filtered 5 kV bank voltage plasma as seen by a 16 channel PDA. The vertical yellow 
lines represent the boundaries of the quiescent period. The red line is the calculated 
emissivity centroid of the plasma. Note that the measured signal drops with increased 
hltering; however, a peaked prohle is still evident in all cases during quiescence. 
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Figure 5.4: Bremsstrahlung emission and line-integrated electron density, (a) shows 
the line-integrated electron number density measured by the interferometer for a 5 kV 
bank voltage plasma at a center (black line) and 1 cm offset (red line), (b) and (c) 
are the square roots of the chord-integrated emission prohles for a bremsstrahlung 
hltered 5 kV bank voltage plasma at centered and approximately 1 cm offset channels 
along the vertical and horizontal axis respectively. The legends indicate the positions 
of the chords. Although the measurements are for different pinches, the shots exhibit 
similar (d) magnetic mode data and quiescent period boundaries, shown as the blue 
and green dashed lines. Both the electron number density and the square root of 
emission show a peaked plasma prohle on axis during quiescence with a lower density 
or emissivity at a 1 cm radius. 
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roughly 500 ns per millimeter and runs approximately 8 cm in length, beginning the 
image during the formation stage allows for the majority of the quiescent period to be 
recorded. Mounted opposite the y-axis viewing PDA at z = 0, the chord-integrated 
emissivity and vertical movement recorded by both diagnostics during quiescence 
is very similar. The centroid of emissivity determined using the theory discussed 
in Section 3.2 also correlates to the camera emission recorded on the him. Slight 
differences in signal amplitude are most likely the result of image saturation and him 
sensitivity of the Imacon photos. Variations in the ampliher calibrations and scaling 
factors also result in diherences between the prohles. Figure 5.5 shows similarities 
between the PDA centroid and Imacon emission prohle for an 8 kV bank voltage 
plasma during a short quiescent period between 29-36 ys. 

Also shown in Figure 5.6 is a 5 kV bank voltage shot with a much longer qui¬ 
escent period. Both of the photos and chord-integrated data are unhltered, and the 
Imacon image is ht to a y-axis similar to the one set for the PDA prohle. With these 
estimations, the PDA and Imacon contours show the same y-axis movement during 
formation and amplitude variations during quiescence. 


5.4 Centroid Results and Comparisons 

Analysis of the center of current measurements using the magnetic probes described in 
Section 2.2 and the emissivity centroid discussed in Section 3.2 conhrms that a direct 
relationship between the two exists during the majority of quiescence. Figure 5.7 
shows an unhltered 5 kV bank voltage plasma viewed in the x-axis by 15 channels 
of a PDA and an 8 probe azimuthal magnetic probe array at the z = 0 plane. The 
azimuthal magnetic probes measure both the amplitude and phase angle associated 
with the m = 1 mode. Projecting these data onto the x-axis and multiplying it 
by the wall radius divided by two gives the x-axis center of current. In general, 
both diagnostic systems assign similar locations for emission and current centers. 
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Figure 5.5: (a) Imacon and (b) PDA views of an unfiltered 8 kV bank voltage plasma. 
Both diagnostics view the pinch through the z = 0 plane’s side windows. Note 
that both systems record the vertical pinch movement. The red line on both plots 
represents the PDA calculated centroid position. 
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Figure 5.6: (a) Imacon and (b) PDA views of an unfiltered 5 kV bank voltage plasma 
showing pinch quiescence. Imacon film saturation and PDA emission scaling result 
in the different emissivity amplitudes. The red line in both plots represents the PDA 
calculated centroid position. 
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Figure 5.7: Emission centroid and magnetic center of current for an unfiltered 5 kV 
bank voltage Z-pinch. The red line represents the PDA calculated emissivity centroid. 
The yellow line is the x-axis component of the m = 1 mode data representing the 
center of current measured by the azimuthal magnetic probes. The vertical yellow 
lines delineate the quiescent period boundaries. 
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Variations between the measurements are due problems with the magnetic probe array 
as well as minor calibration inaccuracies in the PDA amplihers. As the magnetic probe 
data cannot fully resolve higher frequency bulk plasma motion, any quick movement of 
the emission center is not accurately followed by the magnetic probe data; however, 
general comparison shows a direct correlation between the emissivity centroid and 


current center. 
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5.5 Abel Inversion Results 

Axisymmetric 1-D profiles determined using an Abel inversion show flat radial pro¬ 
files during pinch formation, and a peaked radial pinch with a 1 cm radius during 
quiescence. The increasing disparities between plasma emission center and centroid 
position also play an important role in prohle construction. As the plasma emission 
center and centroid diverge, differences between the left and right half path-length 
matrices increase in magnitude, causing discontinuities between the center locations 
of each side. Figure 5.8 is an Abel prohle of unhltered 5 kV bank voltage plasma 
determined by a 15 chord of a PDA viewing along the y-axis of the pinch. During 
formation, the normalized chord-integrated prohle remains peaked while the radial 
prohle remains hat. In quiescence, the radial prohle exhibits a peaked contour with 
an approximately 1 cm radius. A discontinuity arises between the centers of each 
side of the radial emissivities when the centroid algorithm cannot accurately follow 
the plasma contour in an oh-center prohle. However, when the emission centroid ac¬ 
curately follows the plasma center, the inverted prohle exhibits continuity across the 
left and right sides. In general, since the pinch remains near the center during quies¬ 
cence, the inversion prohle at this time accurately models the peaked chord-integrated 
contour. 

5.6 Fourier Deconvolution Results 

The 2-D Fourier technique described in Section 3.10 provides a piecewise constant ra¬ 
dial and cosine varying azimuthal plasma prohle of the ZaP Z-pinch using the theory 
from Section 3.10. Data from unhltered, Ha and bremsstrahlung hltered emissions 
all show movement of the pinch prohle oh axis and amplitude variations during qui¬ 
escence. Two 16 channel PDAs mounted at z = 0 collected chord-integrated data as 
previously mentioned in Section 5.1. The MATLAB code Full2d.m takes these data 
and produces a two-dimensional plasma centered array of emissivity values versus 
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Figure 5.8: Abel inverted axisymmetric radial profiles for an unfiltered 5 kV bank volt¬ 
age Z-pincli. (a-d) The black contour line represents the normalized chord-integrated 
signals along the y-axis. The red and blue lines are the inverted radial prohles for the 
left and right half sides. The points represent the discrete data points. The dashed 
yellow lines in the y-axis contour plot show the boundaries of the quiescent period 
while the solid cyan lines show the time of each shown inverted prohle. (a) During 
formation, the radial prohle of the pinch remains hat. (b) In quiescence however, 
it has a peaked radial prohle with a radius of about 1 cm. (c) When the centroid 
and emission prohle do not accurately match; discontinuity arises between the left 
and right sides of the radial prohle. (d) In cases where the centroid and emission do 
match however, the radial prohle remains continuous across the center. 
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time. A separate vector of x and y centers chosen by the TDerror .m and bestf it .m 
functions combines with the array data to provide lab-frame prohles of the Z-pinch 
versus time. 

Figure 5.9 shows the x- and y-axis chord-integrated profiles of an unhltered 5 kV 
bank voltage plasma while Figure 5.10 shows several deconvolved prohles of this pinch 
at various times during quiescence. Peaked and hollow prohles of the pinch have a 
radius of about 1 cm and exhibit varying position and amplitude. During quiescence, 
movement oh axis is minimal and changes in emission amplitude are smooth. In 
contrast, during the post quiescent period the plasma movement and amplitude both 
change wildly. 

Figures 5.11 and 5.12 are of a similar 5 kV bank voltage plasma viewed through 
the Kodak-Wratten 12 hlter to view Ha emissions. This plasma exhibits a hollow 
prohle throughout the quiescent period with peaks along a radius of about 1 cm from 
the center of the prohle. This is indicative of the Ha shell, where cold neutral gas hit 
by the plasma emits the majority of the radiation seen through the hltered lens. As 
expected, the amplitude of this emission is lower than the unhltered case. 

Figures 5.13 and 5.14 are the chord-integrated prohles and 2-D Fourier prohles 
of a 5 kV bank voltage plasma viewed through the bremsstrahlung hlters listed in 
Section 5.1. As in the unhltered case, both peaked and hollow prohles are seen in 
the deconvolved radial contours, however the signal amplitudes are much lower. The 
radius of these prohles is also about 1 cm in radius and moves oh the ; 2 -axis during 
quiescence. During periods of low chord-integrated signal in the quiescent period, the 
2-D emissivity prohle exhibits a hollow shape, as seen in Figure 5.14(b). When the 
signal increases during a “hash,” emission in the center of the 2-D contour increases, 
resulting in a peaked prohle as shown in Figures 5.14(a),(b),(d). 
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Figure 5.9: x- and ?/-axis contours of an unfiltered 5 kV bank voltage Z-pinch. 16 
chords view the plasma along the y-axis and 15 view it along the x-axis. The dashed 
yellow lines represent the boundaries of the quiescent period. The solid blue lines are 
the times the 2-D data shown in Figure 5.10 were taken. 
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Figure 5.10: 2-D emissivity profiles of the unfiltered 5 kV bank voltage shown in 
Fignre 5.9. (a) and (b) show a partially hollow prohle. The axisymmetric inner shell 
is a prodnct of the constraint pnt on the Fonrier deconvolution analysis discnssed in 
Section 3.11. (c) shows a peaked pinch prohle with a 1 cm radins. (d) is a peaked 
prohle with a 1 cm ohset in the positive x-axis. This occnrs immediately before the 
pinch goes nnstable, as shown by the last vertical bine line in Fignre 5.9. Also note 
that the shells chosen for this deconvolved prohle have larger radii and coarser spacing 
dne to the ohset reducing the nnmber of chords viewing the center of the plasma. 
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Figure 5.11: Ha chord-integrated x- and y-axis contours of a 5 kV bank voltage Z- 
pinch. 16 chords view the plasma along the y-axis and 14 view it along the x-axis. 
The dashed yellow lines represent the boundaries of the quiescent period. The solid 
blue lines are the times 2-D data shown in Figure 5.12 were taken. 
















80 



Figure 5.12: 2-D profiles of Ha emission from the 5 kV profile in Figure 5.11. (a-d) 
Throughout the quiescent period of the pinch, the 2-D radial prohle remains hollow 
with a radius of approximately 1 cm. This is due to collisions between the plasma 
and colder neutral gas at that radius radiating the majority of the Ha radiation. 
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Figure 5.13: Bremsstrahlung chord-integrated x- and y-axis contours of a 5 kV bank 
voltage Z-pinch. 16 chords view the plasma along the y-axis and 15 view it along 
the a:-axis. The dashed yellow lines represent the boundaries of the quiescent period. 
The solid blue lines are the times 2-D data were taken in Figure 5.14. 














82 



Time = 70.254s Time = 75.75)15 

O O 


?3 -2-1 0 1 2 3-3-2-1 0 1 2 3 



X-axis (cm) 

(c) 


X-axis (cm) 

(d) 


Figure 5.14: 2-D profiles of bremsstrahlung emission from the 5 kV profile in Figure 
5.13. (a), (c), and (d) show a peaked prohle with a 1 cm radius during periods of 
bright emission flashes, (b) When the emission is lower during the quiescent period, 
such as indicated by the third blue line in Figure 5.13, the 2-D prohle has a hollow 
shape. Overall emission is also lower compared to the unhltered and Ha case. 
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Chapter 6 

CONCLUSIONS 

Emissivity profiles and tomographic data taken by both PDAs prove valuable in 
understanding the ZaP Flow Z-Pinch. Using the theory and diagnostics discussed in 
previous chapters allows for the comparison of PDA emissivity data to measurements 
made by the Imacon camera, magnetic probes, and interferometer. Verihcation of 
other diagnostic systems by the PDAs increases the conhdence in data interpretation. 
Radial 1-D Abel inverted and 2-D Fourier deconvolved prohles of emissivity are also 
constructed from PDA data, giving a better understanding of the emissivity prohle 
of the ZaP Flow Z-Pinch. Several simple improvements and alternative experimental 
setups are also presented in this chapter. 

6.1 Summary 

Experimental results of chord integrated data show that the emissivity prohle detected 
by the PDAs is peaked during quiescence and matches the visual prohle measured 
by the Imacon streak photos. The PDAs accurately measure amplitude variations 
and movement of the pinch along the viewing axis. !!„ and bremsstrahlung hlters 
signihcantly decrease the detected signal; however, the data still show features of these 
types of plasma emissions during quiescence and are therefore useful for analysis. 

Light hltered for bremsstrahlung radiation for 5 kV bank voltage pinches show 
a peaked prohle on axis similar to interferometry measurements of line integrated 
electron density of 5 kV shots with comparable normalized mode data. This provides 
a general validation that bremsstrahlung radiation is proportional to the square root 
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of the electron number density. 

The centroid of emissivity determined from the centroid.malgorithm also relates 
to the center of current determined from the azimuthal magnetic probes. In general, 
the emission and the centroid appear to overlay during quiescence. As the emission 
moves further off the ; 2 -axis however, the centroid data begins to diverge from the 
viewed peaked emission. This validates the truncation of these points to reduce 
errors in the Abel and Fourier algorithms. 

The Abel inversion algorithm Abel .m correctly deconvolves properly centered ax- 
isymmetric synthetic and experimental radial prohles. Results from 5 kV bank voltage 
plasma show that the Z-pinch maintains a 1 cm peaked and symmetric prohle during 
quiescence. When the center of the inversion prohle signihcantly differs from the cen¬ 
ter of the chord integrated contour, discontinuities arise between the center channels 
of the left and right half sides. 

2-D Fourier deconvolved data determined using the Full2d.mprogram accurately 
measures experimental and synthetic movement off the ; 2 -axis, amplitude variations 
in emission, and both hollow and peaked plasma contours. Experimental 5 kV bank 
voltage plasmas exhibit a radius of about 1 cm during quiescence regardless if they 
are hollow or peaked. The TDerror.m and bestfit.m algorithms accurately offset 
the shell center and change shell radii to produce the best least squares ht of the 
given data. Unhltered data have the largest amplitude and a combination of hollow 
and peaked prohles during quiescence. Ha emissions viewed through a hlter show a 
hollow prohle with a radius of 1 cm, corresponding to the estimated Ha shell size of 
the Z-pinch. Light hltered for bremsstrahlung shows hollow and peaked prohles with 
a lower signal. During times of greater intensity, the emissivity from this radiation 
appears to be more peaked than at other times in the quiescent period, where it 
exhibits a more hollow shape. 
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6.2 Future Work 

Although the photodiode arrays already provide a great deal of information regarding 
the ZaP Flow Z-Pinch plasma, many simple calibrations and changes to experimental 
setup can improve the quantity and quality of the data collected. Problems with 
amplifier calibrations, vacuum tank windows, cabling, and attenuators continuously 
plagued data collection and skewed data interpretation. Although the current design 
is much more robust than the earlier prototype PDA, the electronic circuits were not 
completely reliable. 

Vibrations from mounting and handling the PDA, as well as those from the ex¬ 
periment itself caused circuit-board connections to come loose, causing changes in the 
amplifier responses. Sealing the trimmer capacitors on the amplifier boards shown in 
Figures 4.8 and 4.10 with an epoxy resin should prohibit further damping changes in 
the PDA response. Once completed, new calibration factors should be determined 
using the process described in Section 4.3. The weight of the 36 twisted pair shielded 
cable also puts a great deal of strain on the DB-25 connectors shown in Figure 4.11. 
Proper strain relief has reduced this force, however resoldering of the connector pins 
is required to reduce loose channel connections that appear occasionally during ex¬ 
perimentation. 

Better calibration of the vacuum tank windows will give a much more accurate 
view of the plasma profile. A constant light source with a flat profile over the entire 
viewing area is required to accurately measure window variations and dispersion. 
Combining several opal glass plates will properly flatten a light source but will also 
require a high intensity light to make up for the reduction in signal through each 
plate. 

Reducing variations between the center of current and centroid of emissivity will 
require new calibration factors for the magnetic probes. This recalibration is currently 
underway on the experiment, and simply replacing the calibration values should re- 
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move many of the discrepancies. In addition, the centroid hnder routine can be 
improved to better handle Z-pinch movement further off the center axis. 

Although this thesis presented comparisons between line-integrated electron den¬ 
sity measured by the interferometer and bremsstrahlung emission viewed by the PDA 
for similar plasmas in Section 5.1, no data currently exist where both diagnostics view 
the same local area of the plasma. Removing the side mounted PDA at the z = 0 
plane to permit the interferometer lasers to measure y-axis density will allow for this 
condition. Although the remaining PDA measures emission along the x-axis, this 
should still give a general relationship between bremsstrahlung emission and density 
for the same pulse. 

Improvements in the 2-D algorithm performance will also decrease inaccuracies 
in the resulting deconvolved data. The 2-D prohle center hnder TDerror. m does not 
take the uncertainty of the position into account. While this method, discussed in 
Section 3.11, is satisfactory for this experimentation, a more accurate prohle centering 
can be attained by htting the prohle according to the variance given in the following 
equation, 

[ 56 ^] = [M^] [SE^] + [SM^] [E^] . ( 6 . 1 ) 

6e‘f. is the variance, or square of the standard deviation of the radial prohle. M is the 
position of while is the uncertainty of that measurement. E and SE"^ are the 
chord integrated signals and signal uncertainty respectively. The current algorithm 
only checks the variance associated with the hrst term on the right side. By including 
the second term, a more accurate ht can be expected. 

Up to this point, all data taken with the PDAs have been in one plane of the 
; 2 -axis. Rotating the PDAs so that the chords are oriented along the ; 2 -axis will allow 
for comparison between chord-integrated emissivity and the wavelength of m = 0 and 
m = 1 perturbations during quiescence. Listings of the chord-integrated codes used 
in data analysis of the PDA are given in Appendix A. 
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Appendix A 

MATLAB CHORD INTEGRATING AND ABEL 
INVERSION CODES 

Appendix A includes the MATLAB code and calling functions dealing with chord- 
integrated and Abel inversion data. Explanations of each code are listed in their ap¬ 
propriate subsection. Codes with movie or plotting routines are not included in these 
listings. Shown below is the function calling sequence. The overall code PDAcode.m 
takes data from the MDSplus storage system and creates array variables of the chord 
integrated signals and their time bases. All the subsequent routines handle data from 
both arrays simultaneously, with the exception of the Abell .mand Abel2 .m functions. 


CALLING SEQUENCE 

• PDAcode.m 

— tankaxis.m 

— magneticmode.m 

— centroid, m 

— centroiderror.m 

— abell.m 


abel2.m 
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A.l PDAcode.m 

y. 

yAdam Madson 
y11/04/04 

y 

yPHOTODIODE ARRAY EMISSIVITY CALCULATIONS PROGRAM: PDAcode.m 

y 

yDescription: This is the overall calling program used for analysis of 
y chord integrated and Abel inverted PDA data. Information is 
y called from the MDS plus tree using the mdsopen.m and mdsvalue.m 
y functions. Commenting out various call functions listed below 
y will allow for thecollection of the data of interest at a faster 

y rate. The dataplots.mfunction plots out all chord integrated, 

y centroid, and windowvalidation plots and stores them as .eps or 
y .jpg files. 

y 

y 

ybegin program 
clear all; close all; 
format short; 

yyy user input yyy 

yLOOPING SHOT NUMBER 
"/oshots on given day 
y for i = 3:31 
y shotnum = 50315000+i 

ySINGLE SHOT NUMBER CODES (COMMENT OUT LOOP) 
shotnum = 50316044 

yProbe Location 
probelocation = 0; 

yPROBE ORIENTATION 

y CHANGE VALUES IN tankaxis.m TO REFLECT PDA ORIENTATION 
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y.y.y. mds data access v/1 

yPOSSIBLE LOOP CREATION TO TRUNCATE MDS VALUE BEFORE PLACING IN ARRAY 

yCODE CHANGED TO ALLOW FOR PROCESSING BOTH PDA DATA SIMULTANEOUSLY 
y. - PDAl SPECIFIC ARRAYS WILL END IN "1" 

7. - PDA2 SPECIFIC ARRAYS WILL END IN "2" 

yoMDS access to "alfven" or "landau" computers 

70 [a,status] = mdsopen('landau.aa.Washington.edu: :zapmain\shotnum); 
[a,status] = mdsopenC'allven.aa.Washington.edu::zapmain',shotnum); 


'/.PDA 1 ARRAY GENERATION 

PDAl = zeros (32,16384) ; '/,16384 from 16kB digitzer setting 

Tms_l = zeros(1,16384); 
nonzerol = [] ; 

'/.Storing MDS data to Matlab array 
for i = 1:32 

'/.PDAl Frequency Smoothed Data 

channell = ['offset_baseline(smoothld’ ... 

’ (smoothld(\pdal_' ,num2str (i) \12),13))’ ]; 

PDAl_tree = mdsvalue(channell); 

j = sum(size(PDAl_tree)) ; 
if j>1000 

PDAKi, :) = PDAl.tree^ ; 

end 

end 

j = 0; 

for i = 1:32 

if sum(abs(PDAl(i,1:500))) ~= 0 

nonzerol(i-j) = i; '/.index of row 

else 

j=j+i; 

end 


end 
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nonzerol '/oindex of PDA 1 channels with collected data 

’/oPDA use check 

if size(nonzerol) == [0 0] 

PDAluse = 0; 

else 

PDAluse = 1; 

end 

if PDAluse == 0 

XI = 0; timel=0; Tms_l=0; PDA1=0; nil=0; nl=0; 
Max_Signall=0; 11=0; Max_Signal_valuel=0; 

else 


XI = . 15748.*nonzerol; yoChannel position (cm) for PDAl 

timel = [Mim_of (\pdal_',num2str(nonzerol(1)) ’)’]; 

Tms_l = mdsvalue(timel)'*le6; /(time in microseconds 

'/oSlGNAL OUTPUT FOR PDAl 

PDAl = -PDAl(nonzerol,:); yreturns positive signals 

"/oindex of rows and columns 

"/o row: channel, column: time step 

[ml,nl] = size(PDAl); 

[Max_Signall,II] = max(PDAl); "/oil returns max signal ch. 

Max_Signal_valuel = max(Max_Signall); 

end 


"/oPDA 2 ARRAY GENERATION 

PDA2 = zeros (32, 16384); "/ol6384 from 16kB digitzer setting 

Tms_2 = zeros (1,16384); 
nonzero2 = [] ; 

"/Storing MDS data to Matlab array 
for i = 1:32 

"/oPDA2 Frequency Smoothed Data 
channel2 = ['offset_baseline(smoothld’ 

’ (smoothld(\pda2_',num2str(i) ’,12),13))’ ]; 
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PDA2_tree = mdsvalue(channel2); 

j = sum(size(PDA2_tree)); 
if j>1000 

PDA2(i, :) = PDA2_tree^; 

end 

end 

j = 0; 

for i = 1:32 

if sum(abs(PDA2(i,1:500))) ~= 0 

nonzero2(i-j) = i; "/oindex of row 

else 

end 

end 

nonzero2 7oindex of PDA 2 channels with collected data 

'/oPDA use check 
if size(nonzero2) == [0 0] 

PDA2use = 0; 

else 

PDA2use = 1; 

end 

if PDA2use == 0 

X2 = 0; time2=0; Tms_2=0; PDA2=0; m2=0; n2=0; 

Max_Signal2=0; 12=0; Max_Signal_value2=0; 

else 

X2 = . 15748.*nonzero2; ’/oChannel position (cm) for PDA2 

time2 = [Mim_of (\pda2_',num2str(nonzero2(l)) ’)’]; 

Tms_2 = mdsvalue(time2)'*le6; 7time in microseconds 

'/oSlGNAL OUTPUT FOR PDA2 

PDA2 = -PDA2(nonzero2, :) ; "/oreturns positive signals 

"/olndex of rows and columns 
% row: channel, column: time step 
[m2,n2] = size(PDA2); 
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[Max_Signal2,12] = max(PDA2); °/oI2 returns max signal ch. 

Max_Signal_value2 = max(Max_Signal2); 

end 

'/oWORKSPACE CLEANING 

clear timel time2 status i j channell channel2 a PDAl_tree ... 
PDA2_tree II 12 

y.y.y. tank centered axis v/1 

'/oNOTE: CHANGE VARIABLES ^axisl’ ’axis2’ TO CORRESPOND TO 
% PDA ORIENTATION 

[axisl,axis2] = tankaxis(nonzerol,nonzero2); 

"/oSimple change from array axis to tank axis channel locations 
XI = axisl; 

X2 = axis2; 

“/oWORKSPACE CLEANING 
clear axisl axis2 


y.y.y. magnetic mode data y„y„y„ 

[magtime,mode,model,Cx,Cy,cutoff] = ... 
magneticmode(shotnum,probelocation); 

y.y.y. electron number density comparison y.y.y. 

y.Use on Brem Filter shots and compare to similar pO interferometer 
y. shots E proportional ne''2 

y. interf erometer(PDAl ,Tms_l,XI,nonzerol, shotnum,mode,magtime, cutoff); 


y.y.y. centroid data y.y.y. 

“/.MATRIX AND VECTOR CENTROID METHOD 
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[Centroidl,Centroid2,Tl,T2,PDAlc,PDA2c,dXl,dX2,Centroidlr,... 
Centroid2r,dXlc,dX2c,Centroidlp,Centroid2p,Tlp,T2p] = ... 

centroid(Max_Signall,Max_Signal2,XI,X2,PDA1,PDA2,... 
ml,m2,nl,n2,Tms_l,Tms_2,PDAluse,PDA2use); 


y.y.y. error analysis (delta centroid2) T/X 

[PDAlt,PDA2t,Centroidlt,Centroid2t,Tlt,T2t,dClt,dC2t] = ... 

centroiderror(ml,m2,nl,n2,PDAlc,PDA2c,X1,X2,dXl,dX2,T1,T2,... 
Centroidl,Centroid2,shotnum,Cx,Cy,magtime,PDAluse,PDA2use,... 
dXlc,dX2c); 


y.y.y. abel inversion m 
ypDAi 

[Lshelll,Rshelll] = abell(XI,nonzerol,Tlt,Centroidlt,PDAlt,... 
magtime,model,Cx,Cy, .. . "/oMode data for subplot 

PDAl,shotnum,Tms_l, . . . "/oContour data for subplot 

PDAluse,ml); yPDA in use 


y.PDA2 

[Lshell2,Rshell2] = abel2(X2,nonzero2,T2t,Centroid2t,PDA2t,... 
magtime,model,Cx,Cy, .. . "/oMode data for subplot 

PDA2,shotnum,Tms_2, . . . "/oContour data for subplot 

PDA2use,m2); yPDA in use 


y. DATA PLOTS yxxx. 

dataplots(XI,X2,Centroidl,Centroid2,Tms_l,Tms_2,T1,T2,PDA1,PDA2,... 

Max_Signal_valuel,Max_Signal_value2,nonzerol,nonzero2,shotnum,... 
PDAluse,PDA2use,Cx,Cy,magtime,Tip,T2p,Tit,T2t,Centroidlp,... 
Centroid2p,Centroidlt,Centroid2t,dXlc,dX2c,... 
Centroidlr,Centroid2r); 


'/oLOOPING SHOT NUMBER ENDING 
y end 
























97 


A.2 tankaxis.m 


y. 

yAdam Madson 

yi/18/05 

y 

yARRAY AXIS VALUES CALL FUNCTION: tankaxis.m 

y 

yUescription: This call function properly sets the chord locations 
y of the two photodiode arrays and references them to the proper 
y X- and y-axis contours. 

y 

y 

yUegin program 

function [axisl,axis2] = tankaxis(nonzerol,nonzero2) 
format short 

if size(nonzerol) == [0 0] 

nzl = 1; al = 0; '/"a" is a multiplier to return zero in X axis 

else 

nzl = nonzerol; al = 1; 

end 

if size(nonzero2) == [0 0] 
nz2 = 1; a2 = 0; 

else 

nz2 = nonzero2; a2 = 1; 

end 


yyy Y-axls array locations (cm) °/°IX 

y NOTE: USE ON SIDE MOUNTED PDA’S WHERE CHANNEL 1 IS ON -Y 
y (floor) SIDE 

yraw = [-2.342, -2.184, -2.027, -1.869, -1.712, -1.554, -1.397,... 

-1.240, -1.082, -0.925, -0.767, -0.610, -0.452, -0.295,... 
-0.137, 0.020, 0.178, 0.335, 0.493, 0.650, 0.808,... 

0.965, 1.123, 1.280, 1.438, 1.595, 1.753, 1.910,... 

2.068, 2.225, 2.383, 2.540]; 
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yl = al*yraw; 
y2 = a2*yraw; 


'/oREORDER OF AXIS VALUES FOR 180 DEGREE PDA ROTATION ON TANK 
% NOTE: USE ON SIDE MOUNTED PDA’S WHERE CHANNEL 1 IS ON +Y 
y. (CEILING) SIDE 

ylflip = fliplr(yl); 
y2flip = fliplr(y2); 


yyy X-axIs array locations (cm) V/X 

'/oNOTE: CURRENTLY USING SAME VALUES FOR X AND Y AXIS VIEWS 
y NOTE: -X DIRECTION IS TOWARDS PUMPING SYSTEM ... therefore 
y *USE fliplr(y) VALUES IF CH. 1 IS ON +X SIDE 

y *USE yraw VALUES IF CH.l IS ON -X SIDE 

yUSER CHANGED PARAMETERS BASED ON THAT DAY’S PDA ORIENTATION AND 
y REMOVES UNUSED CHANNELS 

axisl = yl(nzl); yCurrently PDAl @ PO side ... ch.l (-Y) 

axis2 = y2flip(nz2); yCurrently PDA2 @ PO top ... ch.l (+X) 
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A.3 magnet icmode.m 


y. 

yAdam Madson 
yi/09/05 

y 

y.MAGNETIC MODE CALL FUNCTION: magneticmode.m 

y 

yOescription: This call function takes azimuthal magnetic probe 
y data from experiment and determines the normalized mode data 

y and Fourier component center of current data for the x- and y- 

y axis. Various plotting routines relating to different 
y components of this data are listed at the end of the code. 

y 

y 

ybegin program 

function [magtime,mode,model,Cx,Cy,cutoff] = ... 
magneticmode(shotnum,probelocation) 

format short 

yUSER INPUTS IF RUN AS A PROGRAM AND NOT AS A CALL FUNCTION 
y shotnum = 41103039; 

yProbe Location 
y probelocation = 0; 

"/oAccess data from MDS 

y[a,status] = mdsopenC'landau.aa.washington.edu::zapmain',shotnum); 
[a,status] = mdsopenC'alfven.aa.Washington.edu::zapmain\ shotnum); 

yi6384 data points from 16kB digitzer setting 
time = zeros (1,16384) ; 
ml = zeros (1,16384); 
mO = zeros (1,16384) ; 
s = zeros (1,16384) ; 
c = zeros (1,16384) ; 


yWall radius data 
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r_wall=mdsvalue( ’ \r_wall^); 

'/oMagnetic data for positive z-axis placement 

ml = mdsvalueC [ Am_l_p',num2str(probelocation)])'; 

mO = mdsvalueC[Am_0_p',num2str(probelocation)])'; 

magtime = mdsvalue( [Mim_of (\m_l_p\num2str(probelocation) O'])' 

magtime = le6*magtime; 

'/oEmpirical value cutoff at .2 

cutoff = ones(1,length(magtime)); cutoff = .2*cutoff; 

cutoff 1 = 100*(r_wall/2)*cutoff; 

for i = 1:length(mO) 
if m0(i) > .01 

m0(i) = m0(i); 
else m0(i) = .01; 
end 

end 

'/oNormalized mode data 

% Can plot this for axisymmetric plot of center of current 
mode = ml./mO; 

model = 100*(r_wall/2)*mode; 

'/oSine and Cosine Data from Tree 

s = mdsvalueC['\s_ml_p',num2str(probelocation)])'; 
c = mdsvalueC ['\c_ml_p',num2str(probelocation)])'; 

'/oSine and Cosine data corrected for angle to axial 

% probe array location 

sine = sin(atan2(s,c)-157.5*pi/180); 

cosine = cos(atan2(s,c)-157.5*pi/180); 

'/oDisplacement along x and y planes (cm) 

Cx = 100*(r_wall/2)*mode.*cosine; 

Cy = 100*(r_wall/2)*mode.*sine; 

'/oTank circle 

th = (0:127)/128*2*pi; 

X = 5*cos(th); 
y = 5*sin(th); 
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A.4 centroid, m 


y. 

yAdam Madson 
yi/05/05 

y 

yCENTROID CALCULATION CALL FUNCTION: centroid.m 

y 

yUescription: The centroid function determines the centroid of 
y emissivity of the PDA chord integrated signal using either 
y the vector and matrix method (currently used for error analysis) 
y or the trapezoid rule. Centroid information is found for both 
y PDA arrays (designated by the 1 and 2 after stored variables). 
y This program also truncates non-peaked data points and large 
y offsets of the centroid to any one side. 

y 

y 

y begin function 

function [Centroidl,Centroid2,Tl,T2,PDAlc,PDA2c,dXl,dX2,... 

Centroidlr,Centroid2r,dXlc,dX2c,Centroidlp,Centroid2p,... 
Tlp,T2p]... 

= centroid(Max_SignaIl,Max_SignaI2,XI,X2,PDA1,PDA2,... 
ml,m2,nl,n2,Tms_l,Tms_2,PDAluse,PDA2use); 

yCall function if Trapezoid part is uncommented 
yfunction [Centroidl,Tl,PDAlc,dXl,Centroid_Trap] = ... 
y centroid(shotnum,Max_SignaIl,X1,PDA1,ml,nl,Tms_l) 

format short 


yyy centroid calculations using vector and matrix method yyy 

y Used for Error Analysis 

"/oEmissivity with forced zeros on added outside channels 
Etotl = zeros(ml+2,nl); 

Etotl(1,:) = 0; 

Etotl(end,:) = 0; 








102 


Etotl(2:end-1,:) = PDAl; 

Etot2 = zeros(m2+2,n2); 

Etot2(l, : ) = 0; 

Etot2(end,:) = 0; 

Etot2(2:end-1,:) = PDA2; 

'/odelta XI array 

dXlcent = zeros(1 ,ml-l); ’/odelta XI array of inside channels 

dXl = zeros (1 ,ml+l) ; ’/ototal delta XI array 

for i = l:ml-l 

dXlcent(i) = XI(i+1)-XI(i); 

end 

if PDAluse == 0 
dXl = 1; 

else 

dXl(l) = abs(mean(dXlcent)); 
dXl(end) = abs(mean(dXlcent)); 
dXl(2:end-l) = abs(dXlcent); 

end 

'/odelta X2 array 

dX2cent = zeros(l,m2-l); '/delta X2 array of inside channels 

dX2 = zeros (1 ,m2+l) ; '/total delta X2 array 

for i = l:m2-l 

dX2cent(i) = X2(i+1)-X2(i); 

end 

if PDA2use == 0 
dX2 = 1; 

else 

dX2(l) = abs(mean(dX2cent)); 
dX2(end) = abs(mean(dX2cent)); 
dX2(2:end-l) = abs(dX2cent); 

end 

'/Centered XI position calculations 
Xtotl = zeros(1,ml+2); 
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Xlc = zeros(1,ml+l); 

Xtotl(2:end-l) = XI; 

"/oCentered X2 position calculations 
Xtot2 = zeros(1,m2+2); 

X2c = zeros(1,m2+l); 

Xtot2(2:end-1) = X2; 


yoNOTE: The first and last values of Xtot are dependent on the 
% orientation of the PDA on the machine. Make sure that they are 
% to the left and right of the first and last values of X 
% respectively. This assumes that the PDA straddles the zero axis. 

yPDAl 

if Xl(l) <= 0 

Xtotl(l) = Xl(l)-dXl(l); 

else 

Xtotl(l) = Xl(l)+dXl(l); 

end 

if XI(end) >= 0 

Xtotl(end) = XI(end)+dXl(end); 

else 

Xtotl(end) = XI(end)-dXl(end); 

end 

y.PDA2 

if X2(l) <= 0 

Xtot2(l) = X2(l)-dX2(l); 

else 

Xtot2(l) = X2(l)+dX2(l); 

end 

if X2(end) >= 0 

Xtot2(end) = X2(end)+dX2(end); 

else 

Xtot2(end) = X2(end)-dX2(end); 

end 


yCentered XI array 
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for i = l:ml+l 

Xlc(i) = ( Xtotl(i)+Xtotl(i+l) )/2; 

end 

'/oCentered X2 array 
for i = l:m2+l 

X2c(i) = ( Xtot2(i)+Xtot2(i+l) )/2; 

end 

'/oCentered Emissivity signals 
Elc = zeros(ml+1,nl); 

E2c = zeros(m2+l,n2); 

for i = 

for j = l:nl 

Elc(i,j) = ( EtotKi,j)+Etotl(i+l,j) )/2; 

end 

end 

if PDAluse == 0 
Elc = 1; 

end 

for i = l:m2+l 
for j = l:n2 

E2c(i,j) = ( Etot2(i,j)+Etot2(i+l,j) )/2; 

end 

end 

if PDA2use == 0 
E2c = 1; 

end 

'/oCENTROID CALCULATION 
Centroidlr = zeros(l,nl); 

Centroidlr = (dXl.*Xlc*Elc)./(dXl*Elc); 

Centroid2r = zeros(l,n2); 

Centroid2r = (dX2.*X2c*E2c)./(dX2*E2c); 

'/average signal of outside channels 
if PDAluse == 0 
Endavgl = 0; 
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else 

Endavgl = zeros(l,nl); 

Endavgl = (abs(PDAl(1,:))+abs(PDAl(ml,:)))/2; 

end 

if PDA2use == 0 
Endavg2 = 0; 

else 

Endavg2 = zeros(l,n2); 

Endavg2 = (abs(PDA2(1,:))+abs(PDA2(m2,:)))/2; 

end 

yoTruncate PDAl input values to remove non-peaked data sets 
p = 0; '/odummy variable 

if PDAluse == 0 

j=i; 

else 

j=nl; 

end 

for i = l:j 

if Max_Signall(i) >= 3*Endavgl(i) 

P = [p,i]; 

end 

end 

Centroidlp=Centroidlr(p(2:end)); yreturns peaked centroid data 

Tip = Tms_l(p(2:end)); ytime of peaked data 

PDAlp = PDAl (: ,p(2: end)) ; "/osignals at peaked data points 

yiruncate PDA2 input values to remove non-peaked data sets 
p = 0; '/odummy variable 

if PDA2use == 0 

j=i; 

else 

j=n2; 

end 

for i = l:j 

if Max_Signal2(i) >= 3*Endavg2(i) 

P = [p,i]; 
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end 


end 

Centroid2p=Centroid2r (p(2: end)); "/oreturns peaked centroid data 

T2p = Tms_2(p(2: end)); °/otime of peaked data 

PDA2p = PDA2(: ,p(2:end)) ; "/osignals at peaked data points 

'/oTruncate PDAl input values to remove points where centroid is 
'/o all the way to the left or right of the chords 
p = 0; '/odummy variable 

for i = 1:length(Tlp) 

"/oNOTE: Watch inequality change if XI(1) and XI(end) 

% centered around zero 

if Xl(l) <= 0 

if Centroidlp(i) > XI(1) & Centroidlp(i) < XI(end) 

P = [p,i]; 

end 

elseif Xl(l) >= 0 

if Centroidlp(i) < XI(1) & Centroidlp(i) > XI(end) 

P = [p,i]; 

end 

end 

end 

if PDAluse == 0 

Centroidl =0; T1 = 0; PDAlc =0; dXl = 0; 

else 

Centroidl=Centroidlp(p(2:end)); /(returns correct centroid data 
T1 = Tlp(p(2:end)) ; "/otime of truncated data 

PDAlc = PDAlp(: ,p(2: end)); /(signals at correct data points 

end 

'/oTruncate PDA2 input values to remove points where centroid is 
'/o all the way to the left or right of the chords 
p = 0; '/odummy variable 

for i = 1:length(T2p) 

'/oNOTE: Watch inequality change if X2(l) and X2(end) 

'/o centered around zero 

if X2(l) <= 0 

if Centroid2p(i) > X2(l) & Centroid2p(i) < X2(end) 

P = [p,i]; 


end 
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elseif X2(l) >= 0 

if Centroid2p(i) < X2(l) & Centroid2p(i) > X2(end) 

P = [p,i]; 

end 

end 

end 

if PDA2use == 0 

Centroid2 =0; T2 = 0; PDA2c = 0; dX2 = 0; 

else 

Centroid2=Centroid2p(p(2:end)); "/oreturns correct centroid data 
T2 = T2p(p(2: end)); "/otime of truncated data 

PDA2c = PDA2p(: ,p(2: end)) ; 7oSignals at correct data points 

end 

"/oObtain correctly sized center to center dX for centroiderror .m code 
dXlc = [] ; 
dX2c = [] ; 

if PDAluse == 0 
dXlc = 0; 

else 

for i = 1:length(Xlc)-l 

dXlc(i) = Xlc(i+1)-Xlc(i); 

end 

dXlc = abs(dXlc); 

end 

if PDA2use == 0 
dX2c = 0; 

else 

for i = 1:length(X2c)-l 

dX2c(i) = X2c(i+1)-X2c(i); 

end 

dX2c = abs(dX2c); 

end 


y.y.y. centroid calculations using trapezoid rule y.y.y 

“/.NEED TO CHANGE VARIABLE NAMES TO MATCH NEW DESIGNATIONS 
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% Centroid_Trap = zeros(l,nl); 

% dX_Trap= zeros(l,ml-l); 

% Xl_2 = zeros(1; 

% Z = zeros(1,ml); 
y„ for j = l:nl 

y. Z = PDAl(:,j); 

y for i = 

y dX_Trap(i) = XI(i+1)-XI(i); 

y Xl_2(i) = (Xl(i+l)+Xl(i))/2; 

y end 

y 

y Trap_den = zeros(l,ml-l); 

y Trap_num = zeros(l,ml-l); 

y for i = 2:ml 

y Trap_den(i-1) = .5*dX_Trap(i-l)*(Z(i-l)+Z(i)); 

y Trap_num(i-1) = .5*dX_Trap(i-l)*(Z(i-l)+Z(i))*Xl_2(i-l); 

y end 

y Trap_D = sum(Trap_den); 

y Trap_N = sum(Trap_nuin) ; 

y 

y Centroid_Trap(l,j) = Trap_N/Trap_D; 

y end 
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A.5 centroiderror.m 


y. 

yAdam Madson 

yi/13/05 

y 

yCENTROID STANDARD DEVIATION CALL FUNCTION: centroiderror.m 

y 

"/oDescription: The centroid error call function determines the 
y standard deviation of the truncated centroid points determined 
y in the centroid.m call function and passed to PDAcode.m. It 
y plots errorbars of these deviations at every 50th location along 
y the centroid vector and plots them versus time. 

y 

y 

ybegin program 

function [PDAlt,PDA2t,Centroidlt,Centroid2t,Tlt,T2t,dClt,dC2t] = ... 

centroiderror(ml,m2,nl,n2,PDAlc,PDA2c,X1,X2,dXl,dX2,T1,T2,... 
Centroidl,Centroid2,shotnum,Cx,Cy,magtime,... 

PDAluse,PDA2use,dXlc,dX2c) 

format short 

yyy pda i use yyy 

if PDAluse == 0 

Centroidlt = 0; Tit = 0; PDAlt = 0; dClt = 0; 

else 

"/oDelta Signal vector 

dVlc = zeros(ml-2,length(PDAlc)); 

for i = l:ml-2 

for j = 1:length(PDAlc) 

dVlc(i,j) = PDAlc(i+l,j)-PDAlc(i,j); 

end 

end 

dVl = zeros(ml,length(PDAlc)); 
dVl(l, :) = PDAlcd, :) ; 
dVKend, :) = PDAlc (ml, :) ; 
dVl(2:end-1,:) = dVlc; 
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Clden = dXlc*PDAlc; 

dClsquared = zeros(1,length(PDAlc)); 

for i = 1:length(PDAlc) 

dClsquared(i)= (1/(Clden(i).~2) 

(dXlc.*Xl-Centroidl(i) .*dXlc) .~2*dVl(: ,i) .''2; 

end 

dCl = sqrt(dClsquared); 

"/oTruncate to between 20 and 100 mu*s for plotting purposes 
p = 0; '/odummy variable 

for i = l:length(Tl) 

if Tl(i)>=20 & Tl(i)<=100 
P = [p,i]; 

end 

end 

Centroidlt=Centroidl(p(2:end)); 

Tit = T1 (p(2:end)) ; yoTime values for truncated vector 

PDAlt = PDAlc(:,p(2:end)); ^Signals for truncated matrix 
dClt = dCl(p(2:end)) ; "/oDelta centroid at truncated vector 

"/oTruncation of errorbars to sample every 50th location 
y„ from 20-100 mu*s 

Tlsample = zeros(l,floor(length(Tlt)/50)); 
dClsample = zeros(l,floor(length(dClt)/50)); 

Centroidlsample = zeros(l,floor(length(Centroidlt)/50)); 
for i = 1:1:length(Tlsample) 

Centroidlsample(i) = Centroidlt(i*50); 

Tlsample(i) = Tlt(i*50); 
dClsample(i) = dClt(i*50); 

end 

end 

y.y.y. pda 2 use y.y.y 

if PDA2use == 0 

Centroid2t = 0; T2t = 0; PDA2t = 0; dC2t = 0; 

else 

"/oDelta Signal vector 

dV2c = zeros(m2-2,length(PDA2c)); 

for i = l:m2-2 

for j = 1:length(PDA2c) 
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dV2c(i,j) = PDA2c(i+1,j)-PDA2c(i,j); 

end 

end 

dV2 = zeros(m2,length(PDA2c)); 
dV2(l,:) = PDA2c(l, :) ; 
dV2(end,:) = PDA2c(m2,:); 
dV2(2:end-1,:) = dV2c; 

C2den = dX2c*PDA2c; 

dC2squared = zeros(1,length(PDA2c)); 

for i = 1:length(PDA2c) 

dC2squared(i)= (1/(C2den(i).~2) 

(dX2c.*X2-Centroid2(i) .*dX2c) . ~2*dV2(: ,i) .''2; 

end 

dC2 = sqrt(dC2squared); 

"/oTruncate to between 20 and 100 mu*s for plotting purposes 
p = 0; '/odummy variable 

for i = l:length(T2) 

if T2(i)>=20 & T2(i)<=100 
P = [p,i]; 

end 

end 

Centroid2t=Centroid2(p(2:end)); 

T2t = T2(p(2: end)); "/oTime values for truncated vector 

PDA2t = PDA2c(: ,p(2: end)) ; /(Signals for truncated matrix 
dC2t = dC2(p(2: end)) ; "/oDelta centroid at truncated vector 

"/oTruncation of errorbars to sample every 50th location 
% from 20-100 mu*s 

T2sample = zeros(l,floor(length(T2t)/50)); 
dC2sample = zeros(l,floor(length(dC2t)/50)); 

Centroid2sample = zeros(1,floor(length(Centroid2t)/50)); 
for i = 1:1:length(T2sample) 

Centroid2sample(i) = Centroid2t(i*50); 

T2sample(i) = T2t(i*50); 
dC2sample(i) = dC2t(i*50); 

end 


end 
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A.6 abell.m 

y. 

'/oAdam Mads on 
yi/05/05 

y 

'/.ABEL INVERSION CALL FUNCTION: abell.m 

y 

'/.NOTE: Requires centroid.m and centroiderror.m to run. Also uses 
y magneticmode.m data for some plots. 

y 

'/.Description: The Abel routine determines the radial plasma profile 
'/. from the chord integrated data using the Abel inversion method. 
'/. The abell.m and abel2.m routines are identical except that they 
'/. use PDAl or PDA2 data. Movie and plotting routines are listed 
'/. at the end of this code. 

'/. 

'/. 

'/.begin function 

'/. function [Lshelll .Rshelll] = abell (XI ,nonzerol ,Tlt .Centroidlt, . . . 
'/. PDAlt, . . . '/.Chord data used 

'/. magtime,model,Cx,Cy, .. . '/.Mode data for subplot 

'/. PDAl,shotnum,Tms_l, . . . '/.Contour data for subplot 

'/. PDAluse,ml); '/.PDA in use 

'/.function without magnetic data 

'/. function [Lshelll ,Rshelll] = abell (XI ,nonzerol ,Tlt ,Centroidlt, . . . 
'/. PDAlt, . . . '/.Chord data used 

'/. PDAl,shotnum,Tms_l, . . . '/.Contour data for subplot 

'/. PDAluse) ; '/.PDA in use 

format short 

'/.BREAK IF PDA NOT IN USE 
for alpha = 1 

if PDAluse == 0 

Lshelll=0; Rshelll=0; 
break 


else 
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"/ovector and array creation 

n = length(Tlt); 

leftindex = [] ; 

right index = [] ; 

left = [] ; 

right = [] ; 

Lvect = [] ; 

Rvect = [] ; 
dleft = [] ; 
dright = [] ; 

Lshelll = zeros(length(nonzerol),n); 

Rshelll = zeros(length(nonzerol),n); 

"/oAscending order check 
if Xl(l) < Xl(2) 

XI = XI; 

PDAlt = PDAlt; 
elseif Xl(l) > Xl(2) 

XI = fliplr(Xl); 

PDAlt = flipud(PDAlt); 

end 

’/ocreate length vector to define shells by distance from centroid 
for i = l:n 

signal = PDAlt(:,i); 

"/oleft side 

leftindex = find(Xl <= Centroidlt(i)); 
left = XI(leftindex); 

Lvect = abs (left-Centroidlt (i)); "/oChord distances from centroid 
for p = 1:length(Lvect)-l 

dleft(p) = abs(Lvect(p+1)-Lvect(p)); 

end 

dleftavg = mean(dleft); 

"/oShell radii from centroid and 1st chord distance from cent, 
leftvect = [Lvect(l)+dleftavg, Lvect]; 

leftsignal = signal (leftindex); /(signal at each chord 

LI = zeros(length(leftvect)-l); 

L2 = zeros(length(leftvect)-l); 
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for j = 1:lengthdeftvect)-1 

for k = 1:length(leftvect)-l 

if leftvect (k)''2-leftvect (j+1) ~2 > 0 

Ll(j,k) = sqrt(leftvect(k)~2-leftvect(j+1)~2); 

else 

Ll(j,k) = 0; 

end 

end 

for k = 1:length(leftvect)-l 

if lef tvect (k)''2-lef tvect (j) ~2 > 0 

L2(j,k) = sqrt (lef tvect (k+1) ~2-leftvect (j + 1) "'2) ; 

else 

L2(j,k) = 0; 

end 

end 

L = 2*(L1-L2); ’/oleft side length matrix 

end 

LS = L\leftsignal; 

"/oLeft side shell emissivities at given time steps. Row order 
% from outside to nearest to centroid 

Lshelll(1:length(LS),i) = LS; 

"/oright side 

rightindex = find(Xl >= Centroidlt(i)); 
right = XI(rightindex); 

Rvect = abs(right-Centroidlt(i));/(chord distances from centroid 
for p = 1:length(Rvect)-1 

dright(p) = abs(Rvect(p+l)-Rvect(p)); 

end 

drightavg = mean(dright); 

"/oShell radii from centroid and 1st chord distance from cent, 
rightvect = [Rvect, Rvect(end)+drightavg]; 

rightsignal = signal (rightindex); /(signal at each chord 

R1 = zeros(length(rightvect)-l); 

R2 = zeros(length(rightvect)-l); 

for j = 1:length(rightvect)-l 

for k = 1:length(rightvect)-l 
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if rightvect (k+1)''2-rightvect (j) ~2 > 0 

Rl(j,k) = sqrt(rightvect(k+1)~2-rightvect(j)~2); 

else 

Rl(j,k) = 0; 

end 

end 

for k = 1:length(rightvect)-l 

if rightvect (k)''2-rightvect (j) ~2 > 0 

R2(j,k) = sqrt (rightvect (k) ~2-rightvect(j)''2) ; 

else 

R2(j,k) = 0; 

end 

end 

R = 2*(R1-R2); ’/oright side length matrix 

end 

RS = R\rightsignal; 

"/oRight side shell emissivities at given time steps. Row order 

% from outside to nearest to centroid 

Rshelll(length(nonzerol)-length(RS)+l:end,i) = RS; 


end 

“/.PDA USE LOOP 
end 


end 
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Appendix B 

MATLAB 2-D FOURIER DECONVOLUTION CODES 

Appendix B includes the MATLAB code and calling functions dealing with two 
dimensional Fourier deconvolution. Explanations of each code are listed in their ap¬ 
propriate subsection. Codes with movie or plotting routines are not included in these 
listings. Shown below is the function calling sequence. The overall code FULL2d.m 
takes chord integrated data and time variables from the PDAcode.m function and 
determines the deconvolved radial prohle. 


CALLING SEQUENCE 

• FULL2d.m 


— TDerror.m 


bestht.m 
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B.l FULL2d.m 


I 

y. 

yAdam Madson 
7.2/07/05 

y 

7.20 FOURIER DECONVOLUTION FUNCTION: FULL2d.m 

7. 

yDescription: FULL2d is the overall algorithm for the 2-D Fourier 
y deconvolution of PDA emission. Chord integrated data or radial 
y synthetic profiles for both PDAs are called at the beginning of 
y this function. Indeces of the two PDAs are matched before 
y generating the path lengths matricies for each time step. At 
y this point, TDerror.m is called to minimize error by finding 

y the chi"2 value. The final offset at this data point is fed 

y to the bestfit.m routine which returns the radial profile at 
y that time. The program then calls the next set of chord data 

y for the next time step. Once this looping is completed, the 

y plasma radial data is stored in a multidimensional array and the 

y offset locations in vectors. The FULL2d.m program then 

y transforms this profile from the plasma to the lab plane and 

y plots the data. 

y 

y.NOTE: This program takes several minutes to run. Re-routing the 
y call functions so that TDerror.m calls bestfit.m for the 
y majority of its calculations instead of going through its own 
y looping routine might improve calculation speed. 

7 . 

7. 


ybegin program 
y clear all; close all; 
format short; 

y SHOT VALUES 

y FEED IN Xl,X2,PDAlt,PDA2t,Tlt,T2t,Centroidlt,Centroid2t,PDAl, PDA2 
y USE FULL PDA VALUES 
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'/o % CLEARING UNNEEDED VARIABLES 

clear Centroidl Centroid2 Centroidlp Centroid2p Centroidlr ... 
Centroid2r Max_Signall Max_Signal2 Max_Signal_valuel ... 
Max_Signal_value2 PDAlc PDA2c Tip T2p T1 T2 cutoff ... 
dClt dC2t dXl dXlc dX2 dX2c mode model 

warning('off'allO "/osuppress rank deficient warnings that 

% slow down the code 


y. %BIT NOISE ERROR 

Eml = sum(PDAl(:,1:100),2)/100; 

Em2 = sum(PDA2(:,1:100),2)/100; 


'/.INDEX MATCHING OF PDAl AND PDA2 DATA 
m = length(Tlt); 
n = length(T2t); 
storel = [] ; 
store2 = [] ; 


if m <= n 

p = 1; '/.dummy variable 

for i = 1:length(Tlt) 

Index = find(Tlt(i) == T2t(l:end)); 
if length(Index) ~= 0 
storel(l,p) = i; 
store2(l,p) = Index; 
p = P+1; 

end 

end 

Tltt = Tlt(storel); 

PDAltt = PDAlt(:,storel); 

Centroidltt = Centroidlt(storel); 

T2tt = T2t(store2); 

PDA2tt = PDA2t(:,store2); 

Centroid2tt = Centroid2t(store2); 
elseif m > n 

p = 1; '/.dummy variable 

for i = 1:length(T2t) 

Index = find(T2t(i) == Tit(1:end)); 
if length(Index) ~= 0 
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storel(l,p) = Index; 
store2(l,p) = i; 
p = P+1; 

end 

end 

Tltt = Tlt(storel); 

PDAltt = PDAlt(:,storel); 
Centroidltt = Centroidlt(storel); 
T2tt = T2t(store2); 

PDA2tt = PDA2t(:,store2); 
Centroid2tt = Centroid2t(store2); 

end 


y.y.y. path length matrix generation y:/:/, 

yX-Axis (Ex) viewed by PDA on top or bottom view ports 
"/oY-Axis (Ey) viewed by PDA on side view ports 

y. 

“/.NOTE: WANT ORIENTATION FOR PDA PATH LENGTHS TO BE ASCENDING 
y. ORDER IN BOTH TOP AND SIDE VIEWS 

y. 

“/.NOTE: WANT ORIENTATION FOR GRID 1 TO CORRESPOND TO CHORDl OF 
y. PDAl AND PDA2 TO BE LOCATED IN -X,+Y QUADRANT. THEREFORE, 
% TOP PDA NEEDS CHORD LOCATIONS TO GO IN ASCENDING ORDER 
% (-X -> +X) WHILE SIDE PDA NEEDS CHORDS TO GO IN DESCENDING 

y. ORDER (+Y -> -Y) 

yPDA LOCATION 
y. PDAl: XI :PO-side 
% PDA2: X2 :PO-vertical 

if XI (1) > XI (2) "/oSide view (ascending check) 

XI = fliplr(Xl); 

PDAltt = flipud(PDAltt); 

end 

if X2(l) > X2(2) ytop view (ascending check) 

X2 = fliplr(X2); 

PDA2tt = flipud(PDA2tt); 


end 
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'/.USER INPUT 

Ex = PDA2tt; '/.Ex = signal from TOP VIEW 

Ax = X2; 

Centx = Centroid2tt; 

Ey = PDAltt; '/.Ey = signal from SIDE VIEW 

Ay = XI; 

Centy = Centroidltt; 


'/.CENTER OF PROFILE FITTING ROUTINE 

er_array = zeros(15,length(Tltt)); 

Xoffset = zeros(l,length(Tltt)); 

Yoffset = zeros(1,length(Tltt)); 
r_array = zeros(15,length(Tltt)); 

for i = 1:length(Tltt) 

'/.clearing looping variables 
clear offsetX offsetY ErrorMesh 

Xguess = Centx(i); '/.Initial Profile Center guess 

Yguess = Centy(i); 

deltaX = 1.0; 

deltaY =1.0; 

'/.USE FOR SYNTHETIC DATA 
'/. deltaX = .5; 

'/. deltaY = .5; 

dX = zeros(1,length(Ax)-1); 
dY = zeros(1,length(Ay)-1); 

for j = 1:length(Ax)-1 

dX(j) = Ax(j+1)-Ax(j) ; 

end 

for j = 1:length(Ay)-1 

dY(j) = Ay(j+1)-Ay(j); 

end 

dXavg = mean(dX); 
dYavg = mean(dY); 
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Chordavg = (dXavg+dYavg)/2; 

deltaX = .25*length(Ax)*dXavg; 
deltaY = .25*length(Ay)*dYavg; 


'/oLENGTH MATRIX GENERATION 
signal1 = Ex(:,i) ; 
signal2 = Ey(:,i) ; 
signal = [signall; signal2]; 
offsetX = zeros(9,l); 
offsetY = zeros(9,l); 

MinError = [] ; 

ERRl = [] ; 

ERR2 = [] ; 

delta = max([deltaX deltaY]); 

j = 1; 

while delta > .05 

'/.MINIMUM ERROR LOOP IN X AND Y (DEFINED BY DELTA) 

[dEm_s,Xguessnew,Yguessnew,deltaXnew,deltaYnew,min_error,... 
offsetl,offset2,Error,errl,err2,Ecal2] = ... 
TDerror(Ax,Ay,Xguess,Yguess,deltaX,deltaY,... 
signal,Eml,Em2); 

MinError(j) = min_error; 

ERRl(j) = errl; 

ERR2(j) = err2; 
for k = 1:9 

offsetX(k,j) = offsetl(k); 
offsetY(k,j) = offset2(k); 

ErrorMesh(k,j) = Error(k); 

end 

Xguess = Xguessnew; 

Yguess = Yguessnew; 

delta(j) = max([deltaX deltaY]); 


deltaX = deltaXnew; 
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deltaY = deltaYnew; 


if j >= 500 
break 

end 

j = j+1; 

end 

finalerror = min(MinError); 

P = min(find(min(min(ErrorMesh)) == ErrorMesh)); 

Xoff = offsetX(P); 

Yoff = offsetY(P); 

'/oRETURN ARRAY OF BEST RADIAL PROFILES 

[erT,M,Ml,M2,M3,r,X,Y,erraw] = bestfit(Ax,Ay,Xoff,Yoff.signal); 

'/oTIME STEP SAVES 
Xoffset(i) = Xoff; 

Yoffset(i) = Yoff; 

er_array(1:length(erT),i) = erT; 

r_array(1:length(r),i) = r; 

end 


'/oCONVERSION TO POLAR AND CARTESIAN ARRAYS 
er_arr = [] ; 

er_arr = [er_array(l,:);zeros(l,Iength(er_array));... 
zeros(1,Iength(er_array));er_array(2:end,:)]; 

[m n] = size(er_arr); 

q = (0:pi/16: 2*pi); /oNumber of Theta points 

maxZ = zeros(1,length(Tltt)); 

'/oCartesian Mesh with Center radii array generation 
X = zeros(m/3+1,length(q),length(Tltt) ); 
y = zeros(m/3+1,length(q),length(Tltt)); 
z = zeros(m/3+1,length(q),length(Tltt)); 
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xn = zeros(m/3+1,length(q),length(Tltt)); 
yn = zeros(m/3+1,length(q),length(Tltt)); 
zn = zeros(m/3+1,length(q),length(Tltt)); 

for k = 1:length(Tltt) 

"/Polar Coordinate Values 
Rho = zeros(m/3,length(q)); 

Theta = zeros(m/3,length(q)); 

Z = zeros(m/3,length(q)); 

p = 1; 

for i = 1:length(r_array(find(r_array(:,k)),k)) 
for j = l:length(q) 

Thetad, j) = q(j) ; 

Rho(i,j) = r_array(i,k); 

Z(i,j) = er_arr(p,k)+er_arr(p+l,k)*cos(q(j))+... 
er_arr(p+2,k)*sin(q(j)) ; 

end 

P = P+3; 

end 

delta_r_ray=... 

(r_array(length(find(r_array(:,k))),k)-r_array(l,k))/... 

(r_array(length(find(r_array(:,k))),k)-l); 

for i = length(r_array(find(r_array(:,k)),k))+l:m/3 
for j = l:length(q) 

Thetad, j) = q(j) ; 

Rho(i,j) = r_array(i-l,k)+delta_r_ray; 

end 

P = P+3; 

end 

meanZ(k) = mean(Z(l,:)); 

"/CARTESIAN COORDINATE TRANSFORMATION 
"/ only needed if Theta and Rho are vectors 

"/ [th_bogus ,r_bogus] = meshgrid(Theta,Rho); 

"/ [xc,yc] = pol2cart(th_bogus,r_bogus) ; 


[xc(:,:,k),yc(:,:,k),zc(:,:,k)] = pol2cart(Theta,Rho,Z); 
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"/oAdding position at center 

xl = zeros(1,length(q)); 

yl = zeros(1,length(q)); 

zl = ones(1,length(q))*meanZ(k); 

x(l,:,k) = xl; 

y(l,:,k) = yl; 

z(l,:,k) = zl; 

x(2:end,:,k) = xc(:,:,k); 

y(2:end,:,k) = yc(:,:,k); 

z(2:end,:,k) = zc(:,:,k); 

xn(:,:,k) = x(:,:,k)+Xoffset(k); 
yn(:,:,k) = y(:,:,k)+Yoffset(k); 
zn(:,:,k) = z(:,:,k); 

end 
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B.2 TDerror.m 


y. 

yAdam Madson 
7.2/15/05 

y 

'/.TWO DIMENSIONAL FOURIER DECONVOLUTION ERROR FINDER: TDerror.m 

y 

yDescription: TDerror.m is a calling function of TDerror.m. Given 
y an initial guessing range, it will find the minimum error in 
y terms of the chi''2 values of a given profile to minimize error 

y and find the best shell profile center. Data points are picked 

y from a 9 unit grid using the current guess value and distances 
y between guesses. If the minimum returned value of the error is 
y at the center, the routine lowers the spacing of the guess. If 
y it is to the left or right, it moves the Xguess or Yguess value 
y in that direction. 

y 

function [dEm_s,Xguessnew,Yguessnew,deItaXnew,deItaYnew,min_error,... 
offsetl,offset2,Error,errl,err2,EcaI2] = ... 

TDerror(Ax,Ay,Xguess,Yguess,deltaX,deltaY,signal,Eml,Em2); 

format short 

Error = zeros(3,3); 
errl = [] ; 
err2 = [] ; 
min_error = [] ; 

for g = 1:9 

offsetl = [Xguess-deltaX, 

Xguess-deltaX, 

Xguess-deltaX, 

offset2 = [Yguess+deltaY, 

Yguess, 

Yguess-deltaY, 


Xguess, 

Xguess, 

Xguess, 

Yguess+deltaY, 
Yguess, 

Yguess-deltaY, 


Xguess+deltaX;... 
Xguess+deltaX; . . . 
Xguess+deltaX] ; 

Yguess+deltaY;. . . 
Yguess;... 
Yguess-deltaY]; 


er = [] 
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erT = [] ; 

y.X-AXIS CHORDS AND Y-AXIS CHORDS 

'/.NOTE: THESE CHORD LOCATIONS SHOULD BE USED FOR REGULATING 
y. SHELL SIZES 

X = Ax-offsetl(g); 

Y = Ay-offset2(g); 

y.DELTA-X and DELTA-Y vectors 
dX = zeros(l,Iength(X)-l); 
dY = zeros(l,Iength(Y)-l); 

for j = 1:Iength(X)-l 
dX(j) = X(j+1)-X(j) ; 

end 

for j = 1:Iength(Y)-l 
dY(j) = Y(j+1)-Y(j) ; 

end 

dXavg = mean(dX); 
dYavg = mean(dY); 

Chordavg = (dXavg+dYavg)/2; 

ySHELL CREATION (r) 

yoMinimum set by empirical distance of .01 relative from chords 
Xmin = find(abs(X) > .01); 

Xmin = X(Xmin); 

Ymin = find(abs(Y) > .01); 

Ymin = Y(Ymin); 

minimum = [abs(Xmin) abs(Ymin)]; 
innershell = Chordavg + min(minimum); 

"/oMaximum set by 2nd smallest X Y Chord value greater than . 1 
maximum = sort([abs(Xmin(1)) abs(Xmin(end)) ... 

abs(Ymin(l)) abs(Ymin(end))]); 
outershell = Chordavg + maximum(2); 

r = [innershell:Chordavg:outershell] ; "/oShell radii 

Radii = zeros(1,length(r)+l); 

Radii = [0,r]; 
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T1 = zeros(length(X)(length(r)); 

T2 = zeros(length(X)(length(r)); 

51 = zeros(length(Y)(length(r)); 

52 = zeros(length(Y)(length(r)); 

'/oTOP/BOTTOM VIEW 
for j = l:length(X) 

for k = l:length(r) 

if Radii(k+1)''2-X(j) “2 > 0 

Tl(j(k) = sqrt(Radii(k+1)~2-X(j)~2); 

else 

Tl(j(k) = 0; 

end 

end 

for k = l:length(r) 

if Radii(k)"2-X(j)“2 > 0 

T2(j(k) = sqrt (Radii (k)''2-X(j) ~2); 

else 

T2(j(k) = 0; 

end 

end 

T = 2*(T1-T2); 

end 

'/oSIDE/SIDE VIEW 
for j = l:length(Y) 

for k = l:length(r) 

if Radii(k+1)''2-Y(j)''2 > 0 

Sl(j(k) = sqrt(Radii(k+1)~2-Y(j)~2); 

else 

Sl(j(k) = 0; 

end 

end 

for k = l:length(r) 

if Radii(k)"2-Y(j)“2 > 0 

S2(j(k) = sqrt (Radii (k)''2-Y(j) ~2); 

else 

S2(j(k) = 0; 

end 

end 

S = 2*(S1-S2); 
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end 

'/oMIDPOINT GENERATION 

MIDt_pos = (Tl+T2)/2; "/,¥ value midpoints for +Y quadrants II,I 

MIDs_neg = (Sl+S2)/2; "/oX value midpoints for -X quadrants III,II 

MIDt_neg = -MIDt_pos; %¥ value midpoints for -Y quadrants III,IV 

MIDs_pos = -MIDs_neg; "/oX value midpoints for +X quadrants IV,I 

Xn = find(X < 0) ; ’/oindex of negative X values 

Xp = find(X >= 0) ; ’/oindex of positive X values 

Yn = find(Y < 0) ; "/oindex of negative Y values 

Yp = find(Y >= 0) ; 7oindex of positive Y values 

y.THETA VALUES BY QUADRANT FOR TOP AND SIDE VIEWS 
% These values are represented in column order according to 
% the corresponding values of the X and Y vectors respectively 
% EXAMPLE: Any Theta values corresponding to X(l) would be in 

% quadrant II or III due to ascending order of X vector and 

% would be in the first row of the Theta matrix for those 

% quadrants. The midpoint values used would relate to the 

% quadrant sign (+/-) . 

y. 

"/oTop view quadrant I is upper right corner 

Theta_tl = []; ZTop view Theta values for quadrant I 

Theta_tII = [] ; TiTop view Theta values for quadrant II 

Theta_tIII = []; ZTop view Theta values for quadrant III 

Theta_tIV = [] ; TiTop view Theta values for quadrant IV 

"/oSide view quadrant I is upper left corner due to 90 deg 
% theta rotation 

Theta_sl = []; ySide view Theta values for quadrant I 

Theta_sII = []; ySide view Theta values for quadrant II 

Theta_sIII = []; ySide view Theta values for quadrant III 

Theta_sIV = []; ySide view Theta values for quadrant IV 

for k = l:length(r) 

for j = l:length(Xn) 

Theta_tII(j,k) = atan2(MIDt_pos(j,k),X(j)); 

Theta_tIII(j,k) = atan2(MIDt_neg(j,k),X(j)); 

end 

for j = l:length(Xp) 

Theta_tl(j,k) = atan2(MIDt_pos(length(Xn)+j,k),... 
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X(length(Xn)+j)); 

Theta_tIV(j,k) = atan2(MIDt_neg(length(Xn)+j,k),... 
X(length(Xn)+j)); 

end 

for j = l:length(Yn) 

Theta_sIV(j,k) = atan2(Y(j),MIDs_pos(j,k)); 
Theta_sIII(j,k) = atan2(Y(j),MIDs_neg(j,k)); 

end 

for j = l:length(Yp) 

Theta_sl(j,k) = atan2(Y(length(Yn)+j),... 

MIDs_pos(length(Yn)+j,k)); 

Theta_sII(j,k) = atan2(Y(length(Yn)+j),... 
MIDs_neg(length(Yn)+j,k)); 

end 

end 


"/oHemisphere Theta Matricies 
Theta_tN = [Theta_tII; Theta_tl]; 
Theta_tS = [Theta_tIII; Theta_tIV]; 
Theta_sE = [Theta_sIII; Theta_sII]; 
Theta_sW = [Theta_sIV; Theta_sl]; 


"/oTop view upper half 
"/oTop view lower half 
°/oSide view left half 
°/oSide view right half 


'/oSUBMATRlX GENERATION 
yoRadially dependent submatrix 
Ml = [T;S] ; 

[a,b] = size(Ml); 


yCosine dependent submatrix 

topC = [cos(Theta_tN)+cos(Theta_tS)]; 

sideC = [cos(Theta_sE)+cos(Theta_sW)]; 

fullC = [topC;sideC]; 

for j = 1:a 

for k = l:length(r) 

M2(j,k) = (Ml(j,k)/2).*fullC(j,k); 

end 

end 


"/oSine dependent submatrix 
tops = [sin(Theta_tN)+sin(Theta_tS)] ; 
sides = [sin(Theta_sE)+sin(Theta_sW)]; 
fulls = [tops;sides]; 
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for j = 1: a 

for k = l:length(r) 

M3(j,k) = (Ml(j,k)/2).*fullS(j,k); 

end 

end 

'/oFULL GEOMETRIC MATRIX 
M = zeros(a,3*Iength(r)); 

p = 1; 

for j = l:Iength(r) 

M(:,p) = Ml(:,j); 

M(:,p+l) = M2(:,j); 

M(:,p+2) = M3(:,j); 

P = P+3; 

end 

'/oAX I SYMMETRIC 1ST SHELL 

M(:,2:3) = [] ; "/oforce axisymmetric first shell 

[number_E,number_e]=size(M); 

"/oaxisymmetric, cosine, and sine component indicies 

10 = 2:3:number_e; yoaxisymmetric indicies 

11 = 3:3:number_e; ycosine indicies 

12 = 4:3:number_e; ysine indicies 

10= [1,10]; '/oaxisymmetric for truncated M matrix 

'/oLEAST SQUARES FIT 
erraw = zeros(number_e,1); 

Ecal = zeros(number_E,1); 

erraw = M\signal; 

"/oConstraint on Cosine and Sine components 
eraxi = erraw(IO); 

ervar = 2*sqrt(erraw(Il).~2+erraw(I2).~2); 
for j = l:Iength(I2) 

if eraxi(j+l) < ervar(j) 
eraxi(j+l) = ervar(j); 

end 

er(I0(j+l)) = eraxi(j+l); 


end 
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er(1) = eraxi(1); 
er(Il) = erraw(Il); 
er(I2) = erraw(I2); 
erT = er'; 

"/oError calculation 

Ecal = M*erraw; "/oCalculated Signal 

Ecal2 = M*erT; "/oCalculated Signal with Constraint 

yoCHI"2 Error Analysis 
"/oSimulated Data 
dEm_s = [Eml;Em2].“2; 

errl(g) = sum((Ecal-signal).~2./dEm_s); 
err2(g) = sum( (Ecal2-signal) . ~2./dEni_s) ; 

% '/oEQUAL WEIGHTING 

Error(g) = (sum((Ecal-signal).~2./dEm_s)+... 

sum((Ecal2-signal) . ''2./dEm_s))/2; 

end 

yError Minimization 

index_min = find(abs(Error) == min(min(abs(Error)))); 

"/oReduce Delta Variation 
switch index_min(l) 

yShift Xguess and Yguess 
case 1 

Xguessnew = (Xguess-deltaX)/2; 

Yguessnew = (Yguess+deltaY)/2; 
deltaXnew = deltaX*.9; 
deltaYnew = deltaY*.9; 
case 2 

Xguessnew = (Xguess-deltaX)/2; 

Yguessnew = Yguess; 
deltaXnew = deltaX*.9; 
deltaYnew = deltaY*.9; 
case 3 

Xguessnew = (Xguess-deltaX)/2; 

Yguessnew = (Yguess-deltaY)/2; 
deltaXnew = deltaX*.9; 
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deltaYnew = deltaY*.9; 
case 4 

Xguessnew = Xguess; 

Yguessnew = (Yguess+deltaY)/2; 
deltaXnew = deltaX*.9; 
deltaYnew = deltaY*.9; 

'/oCentered DeltaX and DeltaY reduction 
case 5 

Xguessnew = offsetl(5); 

Yguessnew = offset2(5); 
deltaXnew = deltaX*.7; 
deltaYnew = deltaY*.7; 

'/oShift Xguess and Yguess 
case 6 

Xguessnew = Xguess; 

Yguessnew = (Yguess-deltaY)/2; 
deltaXnew = deltaX*.9; 
deltaYnew = deltaY*.9; 
case 7 

Xguessnew = (Xguess+deltaX)/2; 
Yguessnew = (Yguess+deltaY)/2; 
deltaXnew = deltaX*.9; 
deltaYnew = deltaY*.9; 
case 8 

Xguessnew = (Xguess+deltaX)/2; 
Yguessnew = Yguess; 
deltaXnew = deltaX*.9; 
deltaYnew = deltaY*.9; 
case 9 

Xguessnew = (Xguess+deltaX)/2; 
Yguessnew = (Yguess-deltaY)/2; 
deltaXnew = deltaX*.9; 
deltaYnew = deltaY*.9; 
otherwise 

disp('ERRORS 

end 

min_error = min(min(abs(Error))); 
errl = min(errl); 
err2 = min(errl); 
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B.3 bestRt.m 


y. 

yAdam Madson 
7.2/15/05 

y 

'/.TWO DIMENSIONAL FOURIER DECONVOLUTION ERROR FINDER: TDerror.m 

y 

y.Description: Bestfit.m is a calling function of FULL2d.m. It takes 
y the Xoff and Yoff location chosen by the TDerror.m function as 

y the best point to deconvolve the profile. It then creates a 

y path length matrix and performs a least squares fit on the chord 
y integrated data to produce a final radial profile in the plasma 
y reference frame. This data is then fed back into the FULL2d.m 
y code. 

y 

function [erT,M,Ml,M2,M3,r,X,Y,erraw] = ... 
bestfit(Ax,Ay,Xoff,Yoff,signal); 

format short 

Xcenter = Xoff; 

Ycenter = Yoff; 

yX-AXIS CHORDS AND Y-AXIS CHORDS 

'/.NOTE: THESE CHORD LOCATIONS SHOULD BE USED FOR REGULATING 
y SHELL SIZES 

X = Ax-Xcenter; 

Y = Ay-Ycenter; 

y.DELTA-X and DELTA-Y vectors 
dX = zeros(l,Iength(X)-l); 
dY = zeros(l,Iength(Y)-l); 

for j = 1:Iength(X)-l 
dX(j) = X(j+1)-X(j) ; 

end 

for j = 1:Iength(Y)-l 
dY(j) = Y(j+1)-Y(j) ; 
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end 

dXavg = mean(dX); 
dYavg = mean(dY); 

Chordavg = (dXavg+dYavg)/2; 

'/oSHELL CREATION (r) 

'/oMinimum set by an empirical distance of . 1 relative from any chord 
Xmin = find(abs(X) > .01); 

Xmin = X(Xmin); 

Ymin = find(abs(Y) > .01); 

Ymin = Y(Ymin); 

minimum = [abs(Xmin) abs(Ymin)]; 
innershell = Chordavg; + min(minimum); 

'/oMaximum set by 2nd smallest X Y Chord value greater than . 1 
maximum = sort([abs(Xmin(l)) abs(Xmin(end)) ... 

abs(Ymin(l)) abs(Ymin(end))]); 
outershell = Chordavg + maximum(2); 

r = [innershell:Chordavg:outershell] ; yoShell radii 

Radii = zeros(l,length(r)+l); 

Radii = [0,r]; 

T1 = zeros(length(X)(length(r)); 

T2 = zeros(length(X)(length(r)); 

51 = zeros(length(Y)(length(r)); 

52 = zeros(length(Y)(length(r)); 

'/oTOP/BOTTOM VIEW 
for j = l:length(X) 

for k = l:length(r) 

if Radii(k+1)''2-X(j) “2 > 0 

Tl(j(k) = sqrt(Radii(k+1)~2-X(j)~2); 

else 

Tl(j(k) = 0; 

end 

end 

for k = l:length(r) 

if Radii(k)"2-X(j)“2 > 0 

T2(j(k) = sqrt (Radii (k)''2-X(j) ~2); 
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else 

T2(j,k) = 0; 

end 

end 

T = 2*(T1-T2); 

end 

'/oSIDE/SIDE VIEW 
for j = l:length(Y) 

for k = l:length(r) 

if Radii(k+1)''2-Y(j)''2 > 0 

Sl(j,k) = sqrt(Radii(k+1)~2-Y(j)~2); 

else 

Sl(j,k) = 0; 

end 

end 

for k = l:length(r) 

if Radii (k)''2-Y( j) “2 > 0 

S2(j,k) = sqrt (Radii (k)''2-Y(j) ~2); 

else 

S2(j,k) = 0; 

end 

end 

S = 2*(S1-S2); 

end 

'/oMIDPOINT GENERATION 

MIDt_pos = (Tl+T2)/2; ^Y value midpoints for +Y quadrants II,I 

MIDs_neg = (Sl+S2)/2; ’/oX value midpoints for -X quadrants III,II 

MIDt_neg = -MIDt_pos; "/oY value midpoints for -Y quadrants III,IV 

MIDs_pos = -MIDs_neg; "/oX value midpoints for +X quadrants IV,I 

Xn = find(X < 0) ; "/oindex of negative X values 

Xp = find(X >= 0) ; "/oindex of positive X values 

Yn = find(Y < 0) ; "/oindex of negative Y values 

Yp = find(Y >= 0) ; "/oindex of positive Y values 

y.THETA VALUES BY QUADRANT FOR TOP AND SIDE VIEWS 
% These values are represented in column order according to the 
% corresponding values of the X and Y vectors respectively 
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'/o EXAMPLE. Any Theta values corresponding to X(l) would be in 

% quadrant II or III due to ascending order of X vector and would 

% be in the first row of the Theta matrix for those quadrants. 

'/o The midpoint values used would relate to the quadrant sign (+/-) 

y. 

'/oTop view quadrant I is upper right corner 

Theta_tl = []; yiop view Theta values for quadrant I 

Theta_tII = [] ; yiop view Theta values for quadrant II 

Theta_tIII = []; yiop view Theta values for quadrant III 

Theta_tIV = []; yTop view Theta values for quadrant IV 

ySide view quadrant I is upper left corner due to 
y 90 deg theta rotation 

Theta_sl = []; ySide view Theta values for quadrant I 

Theta_sII = []; ySide view Theta values for quadrant II 

Theta_sIII = []; ySide view Theta values for quadrant III 

Theta_sIV = [] ; ySide view Theta values for quadrant IV 

for k = l:length(r) 

for j = l:length(Xn) 

Theta_tII(j,k) = atan2(MIDt_pos(j,k),X(j)); 
y if abs(Theta_tII(j,k)) == pi 

y Theta_tII(j,k) =0; 

y end 

Theta_tIII(j,k) = atan2(MIDt_neg(j,k),X(j)); 
y if abs(Theta_tIII(j,k)) == pi 

y Theta.tllKj ,k) =0; 

y end 

end 

for j = l:length(Xp) 

Theta_tl(j,k) = atan2(MIDt_pos(length(Xn)+j,k),... 
X(length(Xn)+j)); 

Theta_tIV(j,k) = atan2(MIDt_neg(length(Xn)+j,k),... 
X(length(Xn)+j)); 

end 

for j = l:length(Yn) 

Theta_sIV(j,k) = atan2(Y(j),MIDs_pos(j,k)); 

Theta_sIII(j,k) = atan2(Y(j),MIDs_neg(j,k)); 

end 

for j = l:length(Yp) 

Theta_sl(j,k) = atan2(Y(length(Yn)+j),... 
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MIDs_pos(length(Yn)+j,k)); 

Theta_sII(j,k) = atan2(Y(length(Yn)+j),... 
MIDs_neg(length(Yn)+j,k)); 

end 

end 


"/oHemisphere Theta Matricies 
Theta_tN = [Theta_tII; Theta_tl]; 
Theta_tS = [Theta_tIII; Theta_tIV]; 
Theta_sE = [Theta_sIII; Theta_sII]; 
Theta_sW = [Theta_sIV; Theta_sl]; 


°/oTop view upper half 
°/oTop view lower half 
°/oSide view left half 
°/oSide view right half 


“/oSUBMATRIX GENERATION 
yoRadially dependent submatrix 
Ml = [T;S] ; 

[a,b] = size(Ml); 


"/oCosine dependent submatrix 

topC = [cos(Theta_tN)+cos(Theta_tS)] ; 

sideC = [cos(Theta_sE)+cos(Theta_sW)]; 

fullC = [topC;sideC]; 

for j = 1:a 

for k = l:length(r) 

M2(j,k) = (Ml(j,k)/2).*fullC(j,k); 

end 

end 


"/oSine dependent submatrix 

tops = [sin(Theta_tN)+sin(Theta_tS)] ; 

sides = [sin(Theta_sE)+sin(Theta_sW)]; 

fulls = [tops;sides]; 

for j = 1:a 

for k = l:length(r) 

M3(j,k) = (Ml(j,k)/2).*fullS(j,k); 

end 

end 


'/oFULL GEOMETRIC MATRIX 
M = zeros(a,3*length(r)); 


p = 1; 
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for j = l:length(r) 

M(:,p) = Ml(:,j); 

M(:,p+l) = M2(:,j); 

M(:,p+2) = M3(:,j); 

P = P+3; 

end 

'/oAXISYMMETRIC 1ST SHELL 

M(:,2:3) = [] ; "/oforce axisymmetric first shell 

[number_E,number_e]=size(M); 

'/oaxisymmetric, cosine, and sine component indicies 

10 = 2:3:number_e; "/oaxisymmetric indicies 

11 = 3:3:number_e; "/ocosine indicies 

12 = 4:3:number_e; "/sine indicies 

10= [1,10]; '/oaxisymmetric for truncated M matrix 

"/LEAST SQUARES FIT 
erraw = zeros(number_e,1); 

Ecal = zeros(number_E,1); 

erraw = M\signal; 

"/Constraint on Cosine and Sine components 
eraxi = erraw(IO); 

ervar = 2*sqrt(erraw(Il).~2+erraw(I2).~2); 
for j = l:length(I2) 

if eraxi(j+l) < ervar(j) 
eraxi(j+l) = ervar(j); 

end 

er(I0(j+l)) = eraxi(j+l); 

end 

er(1) = eraxi(1); 
er(Il) = erraw(II); 
er(I2) = erraw(I2); 
erT = er'; 



